Actual source code: ex2.c

  1: /*$Id: ex2.c,v 1.25 2001/08/07 21:30:08 bsmith Exp $*/

  3: static char help[] = "Tests MatTranspose(), MatNorm(), MatValid(), and MatAXPY().\n\n";

 5:  #include petscmat.h

  7: #define  TestMatNorm
  8: #define  TestMatTranspose
  9: #define  TestMatNorm_tmat
 10: #define  TestMatAXPY
 11: #undef   TestMatAXPY2

 15: int main(int argc,char **argv)
 16: {
 17:   Mat          mat,tmat = 0;
 18:   int          m = 7,n,i,j,ierr,size,rank,rstart,rend,rect = 0;
 19:   PetscTruth   flg;
 20:   PetscScalar  v, alpha;
 21:   PetscReal    normf,normi,norm1;

 23:   PetscInitialize(&argc,&argv,(char*)0,help);
 24:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
 25:   PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
 26:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 27:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 28:   n = m;
 29:   PetscOptionsHasName(PETSC_NULL,"-rectA",&flg);
 30:   if (flg) {n += 2; rect = 1;}
 31:   PetscOptionsHasName(PETSC_NULL,"-rectB",&flg);
 32:   if (flg) {n -= 2; rect = 1;}

 34:   /* ------- Assemble matrix, test MatValid() --------- */

 36:   MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,n,&mat);
 37:   MatSetFromOptions(mat);
 38:   MatGetOwnershipRange(mat,&rstart,&rend);
 39:   for (i=rstart; i<rend; i++) {
 40:     for (j=0; j<n; j++) {
 41:       v=10*i+j;
 42:       MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
 43:     }
 44:   }
 45:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 46:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

 48:   /* Test whether matrix has been corrupted (just to demonstrate this
 49:      routine) not needed in most application codes. */
 50:   MatValid(mat,(PetscTruth*)&flg);
 51:   if (!flg) SETERRQ(1,"Corrupted matrix.");

 53:   /* ----------------- Test MatNorm()  ----------------- */
 54: #ifdef TestMatNorm
 55:   MatNorm(mat,NORM_FROBENIUS,&normf);
 56:   MatNorm(mat,NORM_1,&norm1);
 57:   MatNorm(mat,NORM_INFINITY,&normi);
 58:   PetscPrintf(PETSC_COMM_WORLD,"original: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",
 59:                      normf,norm1,normi);
 60:   MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
 61: #endif
 62:   /* --------------- Test MatTranspose()  -------------- */
 63: #ifdef TestMatTranspose
 64:   PetscOptionsHasName(PETSC_NULL,"-in_place",&flg);
 65:   if (!rect && flg) {
 66:     MatTranspose(mat,0);   /* in-place transpose */
 67:     tmat = mat; mat = 0;
 68:   } else {      /* out-of-place transpose */
 69:     MatTranspose(mat,&tmat);
 70:   }
 71: #endif
 72:   /* ----------------- Test MatNorm()  ----------------- */
 73: #ifdef TestMatNorm_tmat
 74:   /* Print info about transpose matrix */
 75:   MatNorm(tmat,NORM_FROBENIUS,&normf);
 76:   MatNorm(tmat,NORM_1,&norm1);
 77:   MatNorm(tmat,NORM_INFINITY,&normi);
 78:   PetscPrintf(PETSC_COMM_WORLD,"transpose: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",
 79:                      normf,norm1,normi);
 80:   MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
 81: #endif
 82:   /* ----------------- Test MatAXPY()  ----------------- */
 83: #ifdef TestMatAXPY
 84:   if (mat && !rect) {
 85:     alpha = 1.0;
 86:     PetscOptionsGetScalar(PETSC_NULL,"-alpha",&alpha,PETSC_NULL);
 87:     PetscPrintf(PETSC_COMM_WORLD,"matrix addition:  B = B + alpha * A\n");
 88:     MatAXPY(&alpha,mat,tmat,DIFFERENT_NONZERO_PATTERN);
 89:     MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
 90:   }
 91:   MatDestroy(tmat);
 92: #endif 
 93: #ifdef TestMatAXPY2
 94:     Mat matB;
 95:     alpha = 1.0;
 96:     PetscPrintf(PETSC_COMM_WORLD,"matrix addition:  B = B + alpha * A, SAME_NONZERO_PATTERN\n");
 97:     MatDuplicate(mat,MAT_COPY_VALUES,&matB);
 98:     MatAXPY(&alpha,mat,matB,SAME_NONZERO_PATTERN);
 99:     MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
100:     MatDestroy(matB);

102:     /* get matB that has nonzeros of mat in all even numbers of row and col */
103:     MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,n,&matB);
104:     MatSetFromOptions(matB);
105:     MatGetOwnershipRange(matB,&rstart,&rend);
106:     for (i=rstart; i<rend; i += 2) {
107:       for (j=0; j<n; j += 2) {
108:         v=10*i+j;
109:         MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES);
110:       }
111:     }
112:     MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY);
113:     MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY);
114:     PetscPrintf(PETSC_COMM_WORLD," B: original matrix:\n");
115:     MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
116:     PetscPrintf(PETSC_COMM_WORLD," A(a subset of B):\n");
117:     MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
118:     PetscPrintf(PETSC_COMM_WORLD,"matrix addition:  B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
119:     MatAXPY(&alpha,matB,mat,SUBSET_NONZERO_PATTERN);
120:     MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
121: 
122:     PetscPrintf(PETSC_COMM_WORLD,"matrix addition:  B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
123:     MatAXPY(&alpha,matB,mat,SUBSET_NONZERO_PATTERN);
124:     MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
125: 
126:     MatDestroy(matB);
127: #endif 
128: 

130:   /* Free data structures */
131:   if (mat)  {MatDestroy(mat);}

133:   PetscFinalize();
134:   return 0;
135: }
136: