Actual source code: ex115.c

petsc-3.12.2 2019-11-22
Report Typos and Errors

  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,256,&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:   MatConvert(B,MATIS,MAT_INITIAL_MATRIX,&D);
122:   MatConvert(B,MATIS,MAT_REUSE_MATRIX,&D);
123:   MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
124:   MatConvert(B,MATAIJ,MAT_REUSE_MATRIX,&C);
125:   MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
126:   MatNorm(C,NORM_INFINITY,&err);
127:   if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat AIJ %g",err);
128:   MatDestroy(&C);
129:   MatConvert(D,MATAIJ,MAT_INITIAL_MATRIX,&C);
130:   MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
131:   MatNorm(C,NORM_INFINITY,&err);
132:   if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat IS %g",err);
133:   MatDestroy(&C);
134:   MatDestroy(&D);

136:   /* check MatCreateFromParCSR */
137:   MatHYPREGetParCSR(B,&parcsr);
138:   MatCreateFromParCSR(parcsr,MATAIJ,PETSC_COPY_VALUES,&D);
139:   MatDestroy(&D);
140:   MatCreateFromParCSR(parcsr,MATHYPRE,PETSC_USE_POINTER,&C);

142:   /* check MatMult operations */
143:   MatMultEqual(A,B,4,&flg);
144:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMult B");
145:   MatMultEqual(A,C,4,&flg);
146:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMult C");
147:   MatMultAddEqual(A,B,4,&flg);
148:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultAdd B");
149:   MatMultAddEqual(A,C,4,&flg);
150:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultAdd C");
151:   MatMultTransposeEqual(A,B,4,&flg);
152:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTranspose B");
153:   MatMultTransposeEqual(A,C,4,&flg);
154:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTranspose C");
155:   MatMultTransposeAddEqual(A,B,4,&flg);
156:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTransposeAdd B");
157:   MatMultTransposeAddEqual(A,C,4,&flg);
158:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTransposeAdd C");

160:   /* check PtAP */
161:   if (testptap && M == N) {
162:     Mat pP,hP;

164:     /* PETSc MatPtAP -> output is a MatAIJ
165:        It uses HYPRE functions when -matptap_via hypre is specified at command line */
166:     MatPtAP(A,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pP);
167:     MatPtAP(A,A,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pP);
168:     MatNorm(pP,NORM_INFINITY,&norm);

170:     /* MatPtAP_HYPRE_HYPRE -> output is a MatHYPRE */
171:     MatPtAP(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
172:     MatPtAP(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
173:     MatAXPY(hP,-1.,pP,DIFFERENT_NONZERO_PATTERN);
174:     MatHasOperation(hP,MATOP_NORM,&flg);
175:     if (!flg) { /* TODO add MatNorm_HYPRE */
176:       MatConvert(hP,MATAIJ,MAT_INPLACE_MATRIX,&hP);
177:     }
178:     MatNorm(hP,NORM_INFINITY,&err);
179:     if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatPtAP %g %g",err,norm);
180:     MatDestroy(&hP);

182:     /* MatPtAP_AIJ_HYPRE -> output can be decided at runtime with -matptap_hypre_outtype */
183:     MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
184:     MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
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)A),PETSC_ERR_PLIB,"Error MatPtAP mixed %g %g",err,norm);
192:     MatDestroy(&hP);

194:     MatDestroy(&pP);
195:   }
196:   MatDestroy(&C);
197:   MatDestroy(&B);

199:   /* check MatMatMult */
200:   if (testmatmatmult) {
201:     MatTranspose(A,MAT_INITIAL_MATRIX,&B);
202:     MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&C);
203:     MatConvert(B,MATHYPRE,MAT_INITIAL_MATRIX,&D);

205:     /* PETSc MatMatMult -> output is a MatAIJ
206:        It uses HYPRE functions when -matmatmult_via hypre is specified at command line */
207:     MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pAB);
208:     MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pAB);
209:     MatNorm(pAB,NORM_INFINITY,&norm);

211:     /* MatMatMult_HYPRE_HYPRE -> output is a MatHYPRE */
212:     MatMatMult(C,D,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CD);
213:     MatMatMult(C,D,MAT_REUSE_MATRIX,PETSC_DEFAULT,&CD);
214:     MatAXPY(CD,-1.,pAB,DIFFERENT_NONZERO_PATTERN);
215:     MatHasOperation(CD,MATOP_NORM,&flg);
216:     if (!flg) { /* TODO add MatNorm_HYPRE */
217:       MatConvert(CD,MATAIJ,MAT_INPLACE_MATRIX,&CD);
218:     }
219:     MatNorm(CD,NORM_INFINITY,&err);
220:     if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMatMult %g %g",err,norm);
221:     MatDestroy(&C);
222:     MatDestroy(&D);
223:     MatDestroy(&CD);
224:     MatDestroy(&pAB);

226:     /* When configured with HYPRE, MatMatMatMult is available for the triplet transpose(aij)-aij-aij */
227:     MatCreateTranspose(A,&C);
228:     MatMatMatMult(C,A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CAB);
229:     MatDestroy(&C);
230:     MatTranspose(A,MAT_INITIAL_MATRIX,&C);
231:     MatMatMult(C,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
232:     MatDestroy(&C);
233:     MatMatMult(D,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
234:     MatNorm(C,NORM_INFINITY,&norm);
235:     MatAXPY(C,-1.,CAB,DIFFERENT_NONZERO_PATTERN);
236:     MatNorm(C,NORM_INFINITY,&err);
237:     if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMatMatMult %g %g",err,norm);
238:     MatDestroy(&C);
239:     MatDestroy(&D);
240:     MatDestroy(&CAB);
241:     MatDestroy(&B);
242:   }

244:   /* Check MatView */
245:   MatViewFromOptions(A,NULL,"-view_A");
246:   MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
247:   MatViewFromOptions(B,NULL,"-view_B");

249:   /* Check MatDuplicate/MatCopy */
250:   for (j=0;j<3;j++) {
251:     MatDuplicateOption dop;

253:     dop = MAT_COPY_VALUES;
254:     if (j==1) dop = MAT_DO_NOT_COPY_VALUES;
255:     if (j==2) dop = MAT_SHARE_NONZERO_PATTERN;

257:     for (i=0;i<3;i++) {
258:       MatStructure str;

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

262:       str = DIFFERENT_NONZERO_PATTERN;
263:       if (i==1) str = SAME_NONZERO_PATTERN;
264:       if (i==2) str = SUBSET_NONZERO_PATTERN;

266:       MatDuplicate(A,dop,&C);
267:       MatDuplicate(B,dop,&D);
268:       if (dop != MAT_COPY_VALUES) {
269:         MatCopy(A,C,str);
270:         MatCopy(B,D,str);
271:       }
272:       /* AXPY with AIJ and HYPRE */
273:       MatAXPY(C,-1.0,D,str);
274:       MatNorm(C,NORM_INFINITY,&err);
275:       if (err > PETSC_SMALL) {
276:         MatViewFromOptions(A,NULL,"-view_duplicate_diff");
277:         MatViewFromOptions(B,NULL,"-view_duplicate_diff");
278:         MatViewFromOptions(C,NULL,"-view_duplicate_diff");
279:         SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 1 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
280:       }
281:       /* AXPY with HYPRE and HYPRE */
282:       MatAXPY(D,-1.0,B,str);
283:       if (err > PETSC_SMALL) {
284:         MatViewFromOptions(A,NULL,"-view_duplicate_diff");
285:         MatViewFromOptions(B,NULL,"-view_duplicate_diff");
286:         MatViewFromOptions(D,NULL,"-view_duplicate_diff");
287:         SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 2 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
288:       }
289:       /* Copy from HYPRE to AIJ */
290:       MatCopy(B,C,str);
291:       /* Copy from AIJ to HYPRE */
292:       MatCopy(A,D,str);
293:       /* AXPY with HYPRE and AIJ */
294:       MatAXPY(D,-1.0,C,str);
295:       MatHasOperation(D,MATOP_NORM,&flg);
296:       if (!flg) { /* TODO add MatNorm_HYPRE */
297:         MatConvert(D,MATAIJ,MAT_INPLACE_MATRIX,&D);
298:       }
299:       MatNorm(D,NORM_INFINITY,&err);
300:       if (err > PETSC_SMALL) {
301:         MatViewFromOptions(A,NULL,"-view_duplicate_diff");
302:         MatViewFromOptions(C,NULL,"-view_duplicate_diff");
303:         MatViewFromOptions(D,NULL,"-view_duplicate_diff");
304:         SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 3 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
305:       }
306:       MatDestroy(&C);
307:       MatDestroy(&D);
308:     }
309:   }
310:   MatDestroy(&B);

312:   MatHasCongruentLayouts(A,&flg);
313:   if (flg) {
314:     Vec y,y2;

316:     MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
317:     MatCreateVecs(A,NULL,&y);
318:     MatCreateVecs(B,NULL,&y2);
319:     MatGetDiagonal(A,y);
320:     MatGetDiagonal(B,y2);
321:     VecAXPY(y2,-1.0,y);
322:     VecNorm(y2,NORM_INFINITY,&err);
323:     if (err > PETSC_SMALL) {
324:       VecViewFromOptions(y,NULL,"-view_diagonal_diff");
325:       VecViewFromOptions(y2,NULL,"-view_diagonal_diff");
326:       SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatGetDiagonal %g",err);
327:     }
328:     MatDestroy(&B);
329:     VecDestroy(&y);
330:     VecDestroy(&y2);
331:   }

333:   MatDestroy(&A);

335:   PetscFinalize();
336:   return ierr;
337: }


340: /*TEST

342:    build:
343:       requires: hypre

345:    test:
346:       suffix: 1
347:       args: -N 11 -M 11
348:       output_file: output/ex115_1.out

350:    test:
351:       suffix: 2
352:       nsize: 3
353:       args: -N 13 -M 13 -matmatmult_via hypre
354:       output_file: output/ex115_1.out

356:    test:
357:       suffix: 3
358:       nsize: 4
359:       args: -M 13 -N 7 -matmatmult_via hypre
360:       output_file: output/ex115_1.out

362:    test:
363:       suffix: 4
364:       nsize: 2
365:       args: -M 12 -N 19
366:       output_file: output/ex115_1.out

368:    test:
369:       suffix: 5
370:       nsize: 3
371:       args: -M 13 -N 13 -matptap_via hypre -matptap_hypre_outtype hypre
372:       output_file: output/ex115_1.out

374:    test:
375:       suffix: 6
376:       nsize: 3
377:       args: -M 12 -N 19 -test_offproc
378:       output_file: output/ex115_1.out

380:    test:
381:       suffix: 7
382:       nsize: 3
383:       args: -M 19 -N 12 -test_offproc -view_B ::ascii_info_detail
384:       output_file: output/ex115_7.out

386:    test:
387:       suffix: 8
388:       nsize: 3
389:       args: -M 1 -N 12 -test_offproc
390:       output_file: output/ex115_1.out

392:    test:
393:       suffix: 9
394:       nsize: 3
395:       args: -M 1 -N 2 -test_offproc
396:       output_file: output/ex115_1.out

398: TEST*/