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*/