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