Actual source code: ex53.c


  2: static char help[] = "Tests various routines in MatMPIBAIJ format.\n";

  4: #include <petscmat.h>
  5: #define IMAX 15
  6: int main(int argc,char **args)
  7: {
  8:   Mat               A,B,C,At,Bt;
  9:   PetscViewer       fd;
 10:   char              file[PETSC_MAX_PATH_LEN];
 11:   PetscRandom       rand;
 12:   Vec               xx,yy,s1,s2;
 13:   PetscReal         s1norm,s2norm,rnorm,tol=1.e-10;
 14:   PetscInt          rstart,rend,rows[2],cols[2],m,n,i,j,M,N,ct,row,ncols1,ncols2,bs;
 15:   PetscMPIInt       rank,size;
 16:   PetscErrorCode    0;
 17:   const PetscInt    *cols1,*cols2;
 18:   PetscScalar       vals1[4],vals2[4],v;
 19:   const PetscScalar *v1,*v2;
 20:   PetscBool         flg;

 22:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 23:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 24:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 26:   /* Check out if MatLoad() works */
 27:   PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
 28:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Input file not specified");
 29:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 30:   MatCreate(PETSC_COMM_WORLD,&A);
 31:   MatSetType(A,MATBAIJ);
 32:   MatLoad(A,fd);

 34:   MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);
 35:   PetscViewerDestroy(&fd);

 37:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 38:   PetscRandomSetFromOptions(rand);
 39:   MatGetLocalSize(A,&m,&n);
 40:   VecCreate(PETSC_COMM_WORLD,&xx);
 41:   VecSetSizes(xx,m,PETSC_DECIDE);
 42:   VecSetFromOptions(xx);
 43:   VecDuplicate(xx,&s1);
 44:   VecDuplicate(xx,&s2);
 45:   VecDuplicate(xx,&yy);
 46:   MatGetBlockSize(A,&bs);

 48:   /* Test MatNorm() */
 49:   MatNorm(A,NORM_FROBENIUS,&s1norm);
 50:   MatNorm(B,NORM_FROBENIUS,&s2norm);
 51:   rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
 52:   if (rnorm>tol) {
 53:     PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
 54:   }
 55:   MatNorm(A,NORM_INFINITY,&s1norm);
 56:   MatNorm(B,NORM_INFINITY,&s2norm);
 57:   rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
 58:   if (rnorm>tol) {
 59:     PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
 60:   }
 61:   MatNorm(A,NORM_1,&s1norm);
 62:   MatNorm(B,NORM_1,&s2norm);
 63:   rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
 64:   if (rnorm>tol) {
 65:     PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
 66:   }

 68:   /* Test MatMult() */
 69:   for (i=0; i<IMAX; i++) {
 70:     VecSetRandom(xx,rand);
 71:     MatMult(A,xx,s1);
 72:     MatMult(B,xx,s2);
 73:     VecAXPY(s2,-1.0,s1);
 74:     VecNorm(s2,NORM_2,&rnorm);
 75:     if (rnorm>tol) {
 76:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMult - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
 77:     }
 78:   }

 80:   /* test MatMultAdd() */
 81:   for (i=0; i<IMAX; i++) {
 82:     VecSetRandom(xx,rand);
 83:     VecSetRandom(yy,rand);
 84:     MatMultAdd(A,xx,yy,s1);
 85:     MatMultAdd(B,xx,yy,s2);
 86:     VecAXPY(s2,-1.0,s1);
 87:     VecNorm(s2,NORM_2,&rnorm);
 88:     if (rnorm>tol) {
 89:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultAdd - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
 90:     }
 91:   }

 93:   /* Test MatMultTranspose() */
 94:   for (i=0; i<IMAX; i++) {
 95:     VecSetRandom(xx,rand);
 96:     MatMultTranspose(A,xx,s1);
 97:     MatMultTranspose(B,xx,s2);
 98:     VecNorm(s1,NORM_2,&s1norm);
 99:     VecNorm(s2,NORM_2,&s2norm);
100:     rnorm = s2norm-s1norm;
101:     if (rnorm<-tol || rnorm>tol) {
102:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
103:     }
104:   }
105:   /* Test MatMultTransposeAdd() */
106:   for (i=0; i<IMAX; i++) {
107:     VecSetRandom(xx,rand);
108:     VecSetRandom(yy,rand);
109:     MatMultTransposeAdd(A,xx,yy,s1);
110:     MatMultTransposeAdd(B,xx,yy,s2);
111:     VecNorm(s1,NORM_2,&s1norm);
112:     VecNorm(s2,NORM_2,&s2norm);
113:     rnorm = s2norm-s1norm;
114:     if (rnorm<-tol || rnorm>tol) {
115:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
116:     }
117:   }

119:   /* Check MatGetValues() */
120:   MatGetOwnershipRange(A,&rstart,&rend);
121:   MatGetSize(A,&M,&N);

123:   for (i=0; i<IMAX; i++) {
124:     /* Create random row numbers ad col numbers */
125:     PetscRandomGetValue(rand,&v);
126:     cols[0] = (int)(PetscRealPart(v)*N);
127:     PetscRandomGetValue(rand,&v);
128:     cols[1] = (int)(PetscRealPart(v)*N);
129:     PetscRandomGetValue(rand,&v);
130:     rows[0] = rstart + (int)(PetscRealPart(v)*m);
131:     PetscRandomGetValue(rand,&v);
132:     rows[1] = rstart + (int)(PetscRealPart(v)*m);

134:     MatGetValues(A,2,rows,2,cols,vals1);
135:     MatGetValues(B,2,rows,2,cols,vals2);

137:     for (j=0; j<4; j++) {
138:       if (vals1[j] != vals2[j]) {
139:         PetscPrintf(PETSC_COMM_SELF,"[%d]: Error: MatGetValues rstart = %2d  row = %2d col = %2d val1 = %e val2 = %e bs = %D\n",rank,rstart,rows[j/2],cols[j%2],PetscRealPart(vals1[j]),PetscRealPart(vals2[j]),bs);
140:       }
141:     }
142:   }

144:   /* Test MatGetRow()/ MatRestoreRow() */
145:   for (ct=0; ct<100; ct++) {
146:     PetscRandomGetValue(rand,&v);
147:     row  = rstart + (int)(PetscRealPart(v)*m);
148:     MatGetRow(A,row,&ncols1,&cols1,&v1);
149:     MatGetRow(B,row,&ncols2,&cols2,&v2);

151:     for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
152:       while (cols2[j] != cols1[i]) i++;
153:       if (v1[i] != v2[j]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - vals incorrect.");
154:     }
155:     if (j<ncols2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - cols incorrect");

157:     MatRestoreRow(A,row,&ncols1,&cols1,&v1);
158:     MatRestoreRow(B,row,&ncols2,&cols2,&v2);
159:   }

161:   /* Test MatConvert() */
162:   MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C);

164:   /* See if MatMult Says both are same */
165:   for (i=0; i<IMAX; i++) {
166:     VecSetRandom(xx,rand);
167:     MatMult(A,xx,s1);
168:     MatMult(C,xx,s2);
169:     VecNorm(s1,NORM_2,&s1norm);
170:     VecNorm(s2,NORM_2,&s2norm);
171:     rnorm = s2norm-s1norm;
172:     if (rnorm<-tol || rnorm>tol) {
173:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
174:     }
175:   }
176:   MatDestroy(&C);

178:   /* Test MatTranspose() */
179:   MatTranspose(A,MAT_INITIAL_MATRIX,&At);
180:   MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);
181:   for (i=0; i<IMAX; i++) {
182:     VecSetRandom(xx,rand);
183:     MatMult(At,xx,s1);
184:     MatMult(Bt,xx,s2);
185:     VecNorm(s1,NORM_2,&s1norm);
186:     VecNorm(s2,NORM_2,&s2norm);
187:     rnorm = s2norm-s1norm;
188:     if (rnorm<-tol || rnorm>tol) {
189:       PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
190:     }
191:   }
192:   MatDestroy(&At);
193:   MatDestroy(&Bt);

195:   MatDestroy(&A);
196:   MatDestroy(&B);
197:   VecDestroy(&xx);
198:   VecDestroy(&yy);
199:   VecDestroy(&s1);
200:   VecDestroy(&s2);
201:   PetscRandomDestroy(&rand);
202:   PetscFinalize();
203:   return ierr;
204: }

206: /*TEST

208:    build:
209:       requires: !complex

211:    test:
212:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
213:       nsize: 3
214:       args: -matload_block_size 1 -f ${DATAFILESPATH}/matrices/small

216:    test:
217:       suffix: 2
218:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
219:       nsize: 3
220:       args: -matload_block_size 2 -f ${DATAFILESPATH}/matrices/small
221:       output_file: output/ex53_1.out

223:    test:
224:       suffix: 3
225:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
226:       nsize: 3
227:       args: -matload_block_size 4 -f ${DATAFILESPATH}/matrices/small
228:       output_file: output/ex53_1.out

230:    test:
231:       TODO: Matrix row/column sizes are not compatible with block size
232:       suffix: 4
233:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
234:       nsize: 3
235:       args: -matload_block_size 5 -f ${DATAFILESPATH}/matrices/small
236:       output_file: output/ex53_1.out

238:    test:
239:       suffix: 5
240:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
241:       nsize: 3
242:       args: -matload_block_size 6 -f ${DATAFILESPATH}/matrices/small
243:       output_file: output/ex53_1.out

245:    test:
246:       TODO: Matrix row/column sizes are not compatible with block size
247:       suffix: 6
248:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
249:       nsize: 3
250:       args: -matload_block_size 7 -f ${DATAFILESPATH}/matrices/small
251:       output_file: output/ex53_1.out

253:    test:
254:       TODO: Matrix row/column sizes are not compatible with block size
255:       suffix: 7
256:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
257:       nsize: 3
258:       args: -matload_block_size 8 -f ${DATAFILESPATH}/matrices/small
259:       output_file: output/ex53_1.out

261:    test:
262:       suffix: 8
263:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
264:       nsize: 4
265:       args: -matload_block_size 3 -f ${DATAFILESPATH}/matrices/small
266:       output_file: output/ex53_1.out

268: TEST*/