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:   PetscErrorCode     ierr;
 14:   PetscBool          flg,testptap = PETSC_TRUE,testmatmatmult = PETSC_TRUE;
 15:   PetscReal          norm;
 16:   char               file[256];

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

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

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

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

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

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

101:     /* MatAXPY further exercises MatSetValues_HYPRE */
102:     MatAXPY(B,-1.,A,DIFFERENT_NONZERO_PATTERN);
103:     MatConvert(B,MATMPIAIJ,MAT_INITIAL_MATRIX,&C);
104:     MatNorm(C,NORM_INFINITY,&err);
105:     if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatSetValues %g",err);
106:     MatDestroy(&B);
107:     MatDestroy(&C);
108:   } else {
109:     PetscViewer viewer;

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

118:   /* check conversion routines */
119:   MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
120:   MatConvert(A,MATHYPRE,MAT_REUSE_MATRIX,&B);
121:   MatMultEqual(B,A,4,&flg);
122:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat HYPRE");
123:   MatConvert(B,MATIS,MAT_INITIAL_MATRIX,&D);
124:   MatConvert(B,MATIS,MAT_REUSE_MATRIX,&D);
125:   MatMultEqual(D,A,4,&flg);
126:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat IS");
127:   MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
128:   MatConvert(B,MATAIJ,MAT_REUSE_MATRIX,&C);
129:   MatMultEqual(C,A,4,&flg);
130:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat AIJ");
131:   MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
132:   MatNorm(C,NORM_INFINITY,&err);
133:   if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat AIJ %g",err);
134:   MatDestroy(&C);
135:   MatConvert(D,MATAIJ,MAT_INITIAL_MATRIX,&C);
136:   MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
137:   MatNorm(C,NORM_INFINITY,&err);
138:   if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat IS %g",err);
139:   MatDestroy(&C);
140:   MatDestroy(&D);

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

148:   /* check MatMult operations */
149:   MatMultEqual(A,B,4,&flg);
150:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMult B");
151:   MatMultEqual(A,C,4,&flg);
152:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMult C");
153:   MatMultAddEqual(A,B,4,&flg);
154:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultAdd B");
155:   MatMultAddEqual(A,C,4,&flg);
156:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultAdd C");
157:   MatMultTransposeEqual(A,B,4,&flg);
158:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTranspose B");
159:   MatMultTransposeEqual(A,C,4,&flg);
160:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTranspose C");
161:   MatMultTransposeAddEqual(A,B,4,&flg);
162:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTransposeAdd B");
163:   MatMultTransposeAddEqual(A,C,4,&flg);
164:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTransposeAdd C");

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

170:     /* PETSc MatPtAP -> output is a MatAIJ
171:        It uses HYPRE functions when -matptap_via hypre is specified at command line */
172:     MatPtAP(A,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pP);
173:     MatPtAP(A,A,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pP);
174:     MatNorm(pP,NORM_INFINITY,&norm);
175:     MatPtAPMultEqual(A,A,pP,10,&flg);
176:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP_MatAIJ");

178:     /* MatPtAP_HYPRE_HYPRE -> output is a MatHYPRE */
179:     MatPtAP(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
180:     MatPtAP(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
181:     MatPtAPMultEqual(C,B,hP,10,&flg);
182:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP_HYPRE_HYPRE");

184:     /* Test MatAXPY_Basic() */
185:     MatAXPY(hP,-1.,pP,DIFFERENT_NONZERO_PATTERN);
186:     MatHasOperation(hP,MATOP_NORM,&flg);
187:     if (!flg) { /* TODO add MatNorm_HYPRE */
188:       MatConvert(hP,MATAIJ,MAT_INPLACE_MATRIX,&hP);
189:     }
190:     MatNorm(hP,NORM_INFINITY,&err);
191:     if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)hP),PETSC_ERR_PLIB,"Error MatPtAP %g %g",err,norm);
192:     MatDestroy(&pP);
193:     MatDestroy(&hP);

195:     /* MatPtAP_AIJ_HYPRE -> output can be decided at runtime with -matptap_hypre_outtype */
196:     MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
197:     MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
198:     MatPtAPMultEqual(A,B,hP,10,&flg);
199:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP_AIJ_HYPRE");
200:     MatDestroy(&hP);
201:   }
202:   MatDestroy(&C);
203:   MatDestroy(&B);

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

211:     /* PETSc MatMatMult -> output is a MatAIJ
212:        It uses HYPRE functions when -matmatmult_via hypre is specified at command line */
213:     MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pAB);
214:     MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pAB);
215:     MatNorm(pAB,NORM_INFINITY,&norm);
216:     MatMatMultEqual(A,B,pAB,10,&flg);
217:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatMult_AIJ_AIJ");

219:     /* MatMatMult_HYPRE_HYPRE -> output is a MatHYPRE */
220:     MatMatMult(C,D,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CD);
221:     MatMatMult(C,D,MAT_REUSE_MATRIX,PETSC_DEFAULT,&CD);
222:     MatMatMultEqual(C,D,CD,10,&flg);
223:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatMult_HYPRE_HYPRE");

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

228:     MatHasOperation(CD,MATOP_NORM,&flg);
229:     if (!flg) { /* TODO add MatNorm_HYPRE */
230:       MatConvert(CD,MATAIJ,MAT_INPLACE_MATRIX,&CD);
231:     }
232:     MatNorm(CD,NORM_INFINITY,&err);
233:     if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)CD),PETSC_ERR_PLIB,"Error MatMatMult %g %g",err,norm);

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

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

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

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

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

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

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

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

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

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

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

347:   MatDestroy(&A);

349:   PetscFinalize();
350:   return ierr;
351: }

353: /*TEST

355:    build:
356:       requires: hypre

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

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

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

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

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

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

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

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

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

420: TEST*/