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];

 18:   PetscInitialize(&argc, &args, (char *)0, help);
 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);
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);
123:   MatConvert(B, MATIS, MAT_INITIAL_MATRIX, &D);
124:   MatConvert(B, MATIS, MAT_REUSE_MATRIX, &D);
125:   MatMultEqual(D, A, 4, &flg);
127:   MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &C);
128:   MatConvert(B, MATAIJ, MAT_REUSE_MATRIX, &C);
129:   MatMultEqual(C, A, 4, &flg);
131:   MatAXPY(C, -1., A, SAME_NONZERO_PATTERN);
132:   MatNorm(C, NORM_INFINITY, &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);
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);
151:   MatMultEqual(A, C, 4, &flg);
153:   MatMultAddEqual(A, B, 4, &flg);
155:   MatMultAddEqual(A, C, 4, &flg);
157:   MatMultTransposeEqual(A, B, 4, &flg);
159:   MatMultTransposeEqual(A, C, 4, &flg);
161:   MatMultTransposeAddEqual(A, B, 4, &flg);
163:   MatMultTransposeAddEqual(A, C, 4, &flg);

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);

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);

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);
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);
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);

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);

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);

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);
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: %" PetscInt_FMT " %" PetscInt_FMT "\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:         SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error test 1 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")", 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:         SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error test 2 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")", 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:         SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Error test 3 MatDuplicate/MatCopy %g (%" PetscInt_FMT ",%" PetscInt_FMT ")", 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:       SETERRQ(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 0;
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*/