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:   const PetscInt    *cols1,*cols2;
 17:   PetscScalar       vals1[4],vals2[4],v;
 18:   const PetscScalar *v1,*v2;
 19:   PetscBool         flg;

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

 25:   /* Check out if MatLoad() works */
 26:   PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
 28:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 29:   MatCreate(PETSC_COMM_WORLD,&A);
 30:   MatSetType(A,MATBAIJ);
 31:   MatLoad(A,fd);

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

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

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

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

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

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

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

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

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

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

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

150:     for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
151:       while (cols2[j] != cols1[i]) i++;
153:     }

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

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

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

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

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

205: /*TEST

207:    build:
208:       requires: !complex

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

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

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

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

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

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

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

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

267: TEST*/