Actual source code: ex115.c


  2: static char help[] = "Tests MatHYPRE\n";

  4: #include <petscmathypre.h>

  6: int main(int argc,char **args)
  7: {
  8:   Mat                A,B,C,D;
  9:   Mat                pAB,CD,CAB;
 10:   hypre_ParCSRMatrix *parcsr;
 11:   PetscReal          err;
 12:   PetscInt           i,j,N = 6, M = 6;
 13:   PetscBool          flg,testptap = PETSC_TRUE,testmatmatmult = PETSC_TRUE;
 14:   PetscReal          norm;
 15:   char               file[256];

 17:   PetscInitialize(&argc,&args,(char*)0,help);
 18:   PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
 19: #if defined(PETSC_USE_COMPLEX)
 20:   testptap = PETSC_FALSE;
 21:   testmatmatmult = PETSC_FALSE;
 22:   PetscOptionsInsertString(NULL,"-options_left 0");
 23: #endif
 24:   PetscOptionsGetBool(NULL,NULL,"-ptap",&testptap,NULL);
 25:   PetscOptionsGetBool(NULL,NULL,"-matmatmult",&testmatmatmult,NULL);
 26:   MatCreate(PETSC_COMM_WORLD,&A);
 27:   if (!flg) { /* Create a matrix and test MatSetValues */
 28:     PetscMPIInt size;

 30:     MPI_Comm_size(PETSC_COMM_WORLD,&size);
 31:     PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 32:     PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
 33:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
 34:     MatSetType(A,MATAIJ);
 35:     MatSeqAIJSetPreallocation(A,9,NULL);
 36:     MatMPIAIJSetPreallocation(A,9,NULL,9,NULL);
 37:     MatCreate(PETSC_COMM_WORLD,&B);
 38:     MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);
 39:     MatSetType(B,MATHYPRE);
 40:     if (M == N) {
 41:       MatHYPRESetPreallocation(B,9,NULL,9,NULL);
 42:     } else {
 43:       MatHYPRESetPreallocation(B,6,NULL,6,NULL);
 44:     }
 45:     if (M == N) {
 46:       for (i=0; i<M; i++) {
 47:         PetscInt    cols[] = {0,1,2,3,4,5};
 48:         PetscScalar vals[] = {0,1./size,2./size,3./size,4./size,5./size};
 49:         for (j=i-2; j<i+1; j++) {
 50:           if (j >= N) {
 51:             MatSetValue(A,i,N-1,(1.*j*N+i)/(3.*N*size),ADD_VALUES);
 52:             MatSetValue(B,i,N-1,(1.*j*N+i)/(3.*N*size),ADD_VALUES);
 53:           } else if (i > j) {
 54:             MatSetValue(A,i,PetscMin(j,N-1),(1.*j*N+i)/(2.*N*size),ADD_VALUES);
 55:             MatSetValue(B,i,PetscMin(j,N-1),(1.*j*N+i)/(2.*N*size),ADD_VALUES);
 56:           } else {
 57:             MatSetValue(A,i,PetscMin(j,N-1),-1.-(1.*j*N+i)/(4.*N*size),ADD_VALUES);
 58:             MatSetValue(B,i,PetscMin(j,N-1),-1.-(1.*j*N+i)/(4.*N*size),ADD_VALUES);
 59:           }
 60:         }
 61:         MatSetValues(A,1,&i,PetscMin(6,N),cols,vals,ADD_VALUES);
 62:         MatSetValues(B,1,&i,PetscMin(6,N),cols,vals,ADD_VALUES);
 63:       }
 64:     } else {
 65:       PetscInt  rows[2];
 66:       PetscBool test_offproc = PETSC_FALSE;

 68:       PetscOptionsGetBool(NULL,NULL,"-test_offproc",&test_offproc,NULL);
 69:       if (test_offproc) {
 70:         const PetscInt *ranges;
 71:         PetscMPIInt    rank;

 73:         MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 74:         MatGetOwnershipRanges(A,&ranges);
 75:         rows[0] = ranges[(rank+1)%size];
 76:         rows[1] = ranges[(rank+1)%size + 1];
 77:       } else {
 78:         MatGetOwnershipRange(A,&rows[0],&rows[1]);
 79:       }
 80:       for (i=rows[0];i<rows[1];i++) {
 81:         PetscInt    cols[] = {0,1,2,3,4,5};
 82:         PetscScalar vals[] = {-1,1,-2,2,-3,3};

 84:         MatSetValues(A,1,&i,PetscMin(6,N),cols,vals,INSERT_VALUES);
 85:         MatSetValues(B,1,&i,PetscMin(6,N),cols,vals,INSERT_VALUES);
 86:       }
 87:     }
 88:     /* MAT_FLUSH_ASSEMBLY currently not supported */
 89:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 90:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 91:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 92:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 94: #if defined(PETSC_USE_COMPLEX)
 95:     /* make the matrix imaginary */
 96:     MatScale(A,PETSC_i);
 97:     MatScale(B,PETSC_i);
 98: #endif

100:     /* MatAXPY further exercises MatSetValues_HYPRE */
101:     MatAXPY(B,-1.,A,DIFFERENT_NONZERO_PATTERN);
102:     MatConvert(B,MATMPIAIJ,MAT_INITIAL_MATRIX,&C);
103:     MatNorm(C,NORM_INFINITY,&err);
105:     MatDestroy(&B);
106:     MatDestroy(&C);
107:   } else {
108:     PetscViewer viewer;

110:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&viewer);
111:     MatSetFromOptions(A);
112:     MatLoad(A,viewer);
113:     PetscViewerDestroy(&viewer);
114:     MatGetSize(A,&M,&N);
115:   }

117:   /* check conversion routines */
118:   MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
119:   MatConvert(A,MATHYPRE,MAT_REUSE_MATRIX,&B);
120:   MatMultEqual(B,A,4,&flg);
122:   MatConvert(B,MATIS,MAT_INITIAL_MATRIX,&D);
123:   MatConvert(B,MATIS,MAT_REUSE_MATRIX,&D);
124:   MatMultEqual(D,A,4,&flg);
126:   MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
127:   MatConvert(B,MATAIJ,MAT_REUSE_MATRIX,&C);
128:   MatMultEqual(C,A,4,&flg);
130:   MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
131:   MatNorm(C,NORM_INFINITY,&err);
133:   MatDestroy(&C);
134:   MatConvert(D,MATAIJ,MAT_INITIAL_MATRIX,&C);
135:   MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
136:   MatNorm(C,NORM_INFINITY,&err);
138:   MatDestroy(&C);
139:   MatDestroy(&D);

141:   /* check MatCreateFromParCSR */
142:   MatHYPREGetParCSR(B,&parcsr);
143:   MatCreateFromParCSR(parcsr,MATAIJ,PETSC_COPY_VALUES,&D);
144:   MatDestroy(&D);
145:   MatCreateFromParCSR(parcsr,MATHYPRE,PETSC_USE_POINTER,&C);

147:   /* check MatMult operations */
148:   MatMultEqual(A,B,4,&flg);
150:   MatMultEqual(A,C,4,&flg);
152:   MatMultAddEqual(A,B,4,&flg);
154:   MatMultAddEqual(A,C,4,&flg);
156:   MatMultTransposeEqual(A,B,4,&flg);
158:   MatMultTransposeEqual(A,C,4,&flg);
160:   MatMultTransposeAddEqual(A,B,4,&flg);
162:   MatMultTransposeAddEqual(A,C,4,&flg);

165:   /* check PtAP */
166:   if (testptap && M == N) {
167:     Mat pP,hP;

169:     /* PETSc MatPtAP -> output is a MatAIJ
170:        It uses HYPRE functions when -matptap_via hypre is specified at command line */
171:     MatPtAP(A,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pP);
172:     MatPtAP(A,A,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pP);
173:     MatNorm(pP,NORM_INFINITY,&norm);
174:     MatPtAPMultEqual(A,A,pP,10,&flg);

177:     /* MatPtAP_HYPRE_HYPRE -> output is a MatHYPRE */
178:     MatPtAP(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
179:     MatPtAP(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
180:     MatPtAPMultEqual(C,B,hP,10,&flg);

183:     /* Test MatAXPY_Basic() */
184:     MatAXPY(hP,-1.,pP,DIFFERENT_NONZERO_PATTERN);
185:     MatHasOperation(hP,MATOP_NORM,&flg);
186:     if (!flg) { /* TODO add MatNorm_HYPRE */
187:       MatConvert(hP,MATAIJ,MAT_INPLACE_MATRIX,&hP);
188:     }
189:     MatNorm(hP,NORM_INFINITY,&err);
191:     MatDestroy(&pP);
192:     MatDestroy(&hP);

194:     /* MatPtAP_AIJ_HYPRE -> output can be decided at runtime with -matptap_hypre_outtype */
195:     MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
196:     MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
197:     MatPtAPMultEqual(A,B,hP,10,&flg);
199:     MatDestroy(&hP);
200:   }
201:   MatDestroy(&C);
202:   MatDestroy(&B);

204:   /* check MatMatMult */
205:   if (testmatmatmult) {
206:     MatTranspose(A,MAT_INITIAL_MATRIX,&B);
207:     MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&C);
208:     MatConvert(B,MATHYPRE,MAT_INITIAL_MATRIX,&D);

210:     /* PETSc MatMatMult -> output is a MatAIJ
211:        It uses HYPRE functions when -matmatmult_via hypre is specified at command line */
212:     MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pAB);
213:     MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pAB);
214:     MatNorm(pAB,NORM_INFINITY,&norm);
215:     MatMatMultEqual(A,B,pAB,10,&flg);

218:     /* MatMatMult_HYPRE_HYPRE -> output is a MatHYPRE */
219:     MatMatMult(C,D,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CD);
220:     MatMatMult(C,D,MAT_REUSE_MATRIX,PETSC_DEFAULT,&CD);
221:     MatMatMultEqual(C,D,CD,10,&flg);

224:     /* Test MatAXPY_Basic() */
225:     MatAXPY(CD,-1.,pAB,DIFFERENT_NONZERO_PATTERN);

227:     MatHasOperation(CD,MATOP_NORM,&flg);
228:     if (!flg) { /* TODO add MatNorm_HYPRE */
229:       MatConvert(CD,MATAIJ,MAT_INPLACE_MATRIX,&CD);
230:     }
231:     MatNorm(CD,NORM_INFINITY,&err);

234:     MatDestroy(&C);
235:     MatDestroy(&D);
236:     MatDestroy(&pAB);
237:     MatDestroy(&CD);

239:     /* When configured with HYPRE, MatMatMatMult is available for the triplet transpose(aij)-aij-aij */
240:     MatCreateTranspose(A,&C);
241:     MatMatMatMult(C,A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CAB);
242:     MatDestroy(&C);
243:     MatTranspose(A,MAT_INITIAL_MATRIX,&C);
244:     MatMatMult(C,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
245:     MatDestroy(&C);
246:     MatMatMult(D,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
247:     MatNorm(C,NORM_INFINITY,&norm);
248:     MatAXPY(C,-1.,CAB,DIFFERENT_NONZERO_PATTERN);
249:     MatNorm(C,NORM_INFINITY,&err);
251:     MatDestroy(&C);
252:     MatDestroy(&D);
253:     MatDestroy(&CAB);
254:     MatDestroy(&B);
255:   }

257:   /* Check MatView */
258:   MatViewFromOptions(A,NULL,"-view_A");
259:   MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
260:   MatViewFromOptions(B,NULL,"-view_B");

262:   /* Check MatDuplicate/MatCopy */
263:   for (j=0;j<3;j++) {
264:     MatDuplicateOption dop;

266:     dop = MAT_COPY_VALUES;
267:     if (j==1) dop = MAT_DO_NOT_COPY_VALUES;
268:     if (j==2) dop = MAT_SHARE_NONZERO_PATTERN;

270:     for (i=0;i<3;i++) {
271:       MatStructure str;

273:       PetscPrintf(PETSC_COMM_WORLD,"Dup/Copy tests: %" PetscInt_FMT " %" PetscInt_FMT "\n",j,i);

275:       str = DIFFERENT_NONZERO_PATTERN;
276:       if (i==1) str = SAME_NONZERO_PATTERN;
277:       if (i==2) str = SUBSET_NONZERO_PATTERN;

279:       MatDuplicate(A,dop,&C);
280:       MatDuplicate(B,dop,&D);
281:       if (dop != MAT_COPY_VALUES) {
282:         MatCopy(A,C,str);
283:         MatCopy(B,D,str);
284:       }
285:       /* AXPY with AIJ and HYPRE */
286:       MatAXPY(C,-1.0,D,str);
287:       MatNorm(C,NORM_INFINITY,&err);
288:       if (err > PETSC_SMALL) {
289:         MatViewFromOptions(A,NULL,"-view_duplicate_diff");
290:         MatViewFromOptions(B,NULL,"-view_duplicate_diff");
291:         MatViewFromOptions(C,NULL,"-view_duplicate_diff");
292:         SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 1 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")",err,j,i);
293:       }
294:       /* AXPY with HYPRE and HYPRE */
295:       MatAXPY(D,-1.0,B,str);
296:       if (err > PETSC_SMALL) {
297:         MatViewFromOptions(A,NULL,"-view_duplicate_diff");
298:         MatViewFromOptions(B,NULL,"-view_duplicate_diff");
299:         MatViewFromOptions(D,NULL,"-view_duplicate_diff");
300:         SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 2 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")",err,j,i);
301:       }
302:       /* Copy from HYPRE to AIJ */
303:       MatCopy(B,C,str);
304:       /* Copy from AIJ to HYPRE */
305:       MatCopy(A,D,str);
306:       /* AXPY with HYPRE and AIJ */
307:       MatAXPY(D,-1.0,C,str);
308:       MatHasOperation(D,MATOP_NORM,&flg);
309:       if (!flg) { /* TODO add MatNorm_HYPRE */
310:         MatConvert(D,MATAIJ,MAT_INPLACE_MATRIX,&D);
311:       }
312:       MatNorm(D,NORM_INFINITY,&err);
313:       if (err > PETSC_SMALL) {
314:         MatViewFromOptions(A,NULL,"-view_duplicate_diff");
315:         MatViewFromOptions(C,NULL,"-view_duplicate_diff");
316:         MatViewFromOptions(D,NULL,"-view_duplicate_diff");
317:         SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 3 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")",err,j,i);
318:       }
319:       MatDestroy(&C);
320:       MatDestroy(&D);
321:     }
322:   }
323:   MatDestroy(&B);

325:   MatHasCongruentLayouts(A,&flg);
326:   if (flg) {
327:     Vec y,y2;

329:     MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
330:     MatCreateVecs(A,NULL,&y);
331:     MatCreateVecs(B,NULL,&y2);
332:     MatGetDiagonal(A,y);
333:     MatGetDiagonal(B,y2);
334:     VecAXPY(y2,-1.0,y);
335:     VecNorm(y2,NORM_INFINITY,&err);
336:     if (err > PETSC_SMALL) {
337:       VecViewFromOptions(y,NULL,"-view_diagonal_diff");
338:       VecViewFromOptions(y2,NULL,"-view_diagonal_diff");
339:       SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatGetDiagonal %g",err);
340:     }
341:     MatDestroy(&B);
342:     VecDestroy(&y);
343:     VecDestroy(&y2);
344:   }

346:   MatDestroy(&A);

348:   PetscFinalize();
349:   return 0;
350: }

352: /*TEST

354:    build:
355:       requires: hypre

357:    test:
358:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
359:       suffix: 1
360:       args: -N 11 -M 11
361:       output_file: output/ex115_1.out

363:    test:
364:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
365:       suffix: 2
366:       nsize: 3
367:       args: -N 13 -M 13 -matmatmult_via hypre
368:       output_file: output/ex115_1.out

370:    test:
371:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
372:       suffix: 3
373:       nsize: 4
374:       args: -M 13 -N 7 -matmatmult_via hypre
375:       output_file: output/ex115_1.out

377:    test:
378:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
379:       suffix: 4
380:       nsize: 2
381:       args: -M 12 -N 19
382:       output_file: output/ex115_1.out

384:    test:
385:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
386:       suffix: 5
387:       nsize: 3
388:       args: -M 13 -N 13 -matptap_via hypre -matptap_hypre_outtype hypre
389:       output_file: output/ex115_1.out

391:    test:
392:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
393:       suffix: 6
394:       nsize: 3
395:       args: -M 12 -N 19 -test_offproc
396:       output_file: output/ex115_1.out

398:    test:
399:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
400:       suffix: 7
401:       nsize: 3
402:       args: -M 19 -N 12 -test_offproc -view_B ::ascii_info_detail
403:       output_file: output/ex115_7.out

405:    test:
406:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
407:       suffix: 8
408:       nsize: 3
409:       args: -M 1 -N 12 -test_offproc
410:       output_file: output/ex115_1.out

412:    test:
413:       requires: !defined(PETSC_HAVE_HYPRE_DEVICE)
414:       suffix: 9
415:       nsize: 3
416:       args: -M 1 -N 2 -test_offproc
417:       output_file: output/ex115_1.out

419: TEST*/