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