Actual source code: ex115.c
petsc-3.12.0 2019-09-29
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*/