Actual source code: ex115.c
petsc-3.11.3 2019-06-26
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;
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: MatCreate(PETSC_COMM_WORLD,&A);
21: if (!flg) { /* Create a matrix and test MatSetValues */
22: PetscMPIInt size;
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
26: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
27: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
28: MatSetType(A,MATAIJ);
29: MatSeqAIJSetPreallocation(A,9,NULL);
30: MatMPIAIJSetPreallocation(A,9,NULL,9,NULL);
31: MatCreate(PETSC_COMM_WORLD,&B);
32: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);
33: MatSetType(B,MATHYPRE);
34: if (M == N) {
35: MatHYPRESetPreallocation(B,9,NULL,9,NULL);
36: } else {
37: MatHYPRESetPreallocation(B,6,NULL,6,NULL);
38: }
39: if (M == N) {
40: for (i=0; i<M; i++) {
41: PetscInt cols[] = {0,1,2,3,4,5};
42: PetscScalar vals[] = {0,1./size,2./size,3./size,4./size,5./size};
43: for (j=i-2; j<i+1; j++) {
44: if (j >= N) {
45: MatSetValue(A,i,N-1,(1.*j*N+i)/(3.*N*size),ADD_VALUES);
46: MatSetValue(B,i,N-1,(1.*j*N+i)/(3.*N*size),ADD_VALUES);
47: } else if (i > j) {
48: MatSetValue(A,i,PetscMin(j,N-1),(1.*j*N+i)/(2.*N*size),ADD_VALUES);
49: MatSetValue(B,i,PetscMin(j,N-1),(1.*j*N+i)/(2.*N*size),ADD_VALUES);
50: } else {
51: MatSetValue(A,i,PetscMin(j,N-1),-1.-(1.*j*N+i)/(4.*N*size),ADD_VALUES);
52: MatSetValue(B,i,PetscMin(j,N-1),-1.-(1.*j*N+i)/(4.*N*size),ADD_VALUES);
53: }
54: }
55: MatSetValues(A,1,&i,PetscMin(6,N),cols,vals,ADD_VALUES);
56: MatSetValues(B,1,&i,PetscMin(6,N),cols,vals,ADD_VALUES);
57: }
58: } else {
59: PetscInt rows[2];
60: PetscBool test_offproc = PETSC_FALSE;
62: PetscOptionsGetBool(NULL,NULL,"-test_offproc",&test_offproc,NULL);
63: if (test_offproc) {
64: const PetscInt *ranges;
65: PetscMPIInt rank;
67: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
68: MatGetOwnershipRanges(A,&ranges);
69: rows[0] = ranges[(rank+1)%size];
70: rows[1] = ranges[(rank+1)%size + 1];
71: } else {
72: MatGetOwnershipRange(A,&rows[0],&rows[1]);
73: }
74: for (i=rows[0];i<rows[1];i++) {
75: PetscInt cols[] = {0,1,2,3,4,5};
76: PetscScalar vals[] = {-1,1,-2,2,-3,3};
78: MatSetValues(A,1,&i,PetscMin(6,N),cols,vals,INSERT_VALUES);
79: MatSetValues(B,1,&i,PetscMin(6,N),cols,vals,INSERT_VALUES);
80: }
81: }
82: /* MAT_FLUSH_ASSEMBLY currently not supported */
83: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
85: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
88: /* MatAXPY further exercises MatSetValues_HYPRE */
89: MatAXPY(B,-1.,A,DIFFERENT_NONZERO_PATTERN);
90: MatConvert(B,MATMPIAIJ,MAT_INITIAL_MATRIX,&C);
91: MatNorm(C,NORM_INFINITY,&err);
92: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatSetValues %g",err);
93: MatDestroy(&B);
94: MatDestroy(&C);
95: } else {
96: PetscViewer viewer;
98: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&viewer);
99: MatSetFromOptions(A);
100: MatLoad(A,viewer);
101: PetscViewerDestroy(&viewer);
102: MatGetSize(A,&M,&N);
103: }
105: /* check conversion routines */
106: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
107: MatConvert(A,MATHYPRE,MAT_REUSE_MATRIX,&B);
108: MatConvert(B,MATIS,MAT_INITIAL_MATRIX,&D);
109: MatConvert(B,MATIS,MAT_REUSE_MATRIX,&D);
110: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
111: MatConvert(B,MATAIJ,MAT_REUSE_MATRIX,&C);
112: MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
113: MatNorm(C,NORM_INFINITY,&err);
114: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat AIJ %g",err);
115: MatDestroy(&C);
116: MatConvert(D,MATAIJ,MAT_INITIAL_MATRIX,&C);
117: MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
118: MatNorm(C,NORM_INFINITY,&err);
119: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat IS %g",err);
120: MatDestroy(&C);
121: MatDestroy(&D);
123: /* check MatCreateFromParCSR */
124: MatHYPREGetParCSR(B,&parcsr);
125: MatCreateFromParCSR(parcsr,MATAIJ,PETSC_COPY_VALUES,&D);
126: MatDestroy(&D);
127: MatCreateFromParCSR(parcsr,MATHYPRE,PETSC_USE_POINTER,&C);
129: /* check MatMult operations */
130: MatMultEqual(A,B,4,&flg);
131: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMult B");
132: MatMultEqual(A,C,4,&flg);
133: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMult C");
134: MatMultAddEqual(A,B,4,&flg);
135: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultAdd B");
136: MatMultAddEqual(A,C,4,&flg);
137: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultAdd C");
138: MatMultTransposeEqual(A,B,4,&flg);
139: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTranspose B");
140: MatMultTransposeEqual(A,C,4,&flg);
141: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTranspose C");
142: MatMultTransposeAddEqual(A,B,4,&flg);
143: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTransposeAdd B");
144: MatMultTransposeAddEqual(A,C,4,&flg);
145: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTransposeAdd C");
147: /* check PtAP */
148: if (M == N) {
149: Mat pP,hP;
151: /* PETSc MatPtAP -> output is a MatAIJ
152: It uses HYPRE functions when -matptap_via hypre is specified at command line */
153: MatPtAP(A,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pP);
154: MatPtAP(A,A,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pP);
155: MatNorm(pP,NORM_INFINITY,&norm);
157: /* MatPtAP_HYPRE_HYPRE -> output is a MatHYPRE */
158: MatPtAP(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
159: MatPtAP(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
160: MatAXPY(hP,-1.,pP,DIFFERENT_NONZERO_PATTERN);
161: MatHasOperation(hP,MATOP_NORM,&flg);
162: if (!flg) { /* TODO add MatNorm_HYPRE */
163: MatConvert(hP,MATAIJ,MAT_INPLACE_MATRIX,&hP);
164: }
165: MatNorm(hP,NORM_INFINITY,&err);
166: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatPtAP %g %g",err,norm);
167: MatDestroy(&hP);
169: /* MatPtAP_AIJ_HYPRE -> output can be decided at runtime with -matptap_hypre_outtype */
170: MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
171: MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
172: MatAXPY(hP,-1.,pP,DIFFERENT_NONZERO_PATTERN);
173: MatHasOperation(hP,MATOP_NORM,&flg);
174: if (!flg) { /* TODO add MatNorm_HYPRE */
175: MatConvert(hP,MATAIJ,MAT_INPLACE_MATRIX,&hP);
176: }
177: MatNorm(hP,NORM_INFINITY,&err);
178: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatPtAP mixed %g %g",err,norm);
179: MatDestroy(&hP);
181: MatDestroy(&pP);
182: }
183: MatDestroy(&C);
184: MatDestroy(&B);
186: /* check MatMatMult */
187: MatTranspose(A,MAT_INITIAL_MATRIX,&B);
188: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&C);
189: MatConvert(B,MATHYPRE,MAT_INITIAL_MATRIX,&D);
191: /* PETSc MatMatMult -> output is a MatAIJ
192: It uses HYPRE functions when -matmatmult_via hypre is specified at command line */
193: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pAB);
194: MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pAB);
195: MatNorm(pAB,NORM_INFINITY,&norm);
197: /* MatMatMult_HYPRE_HYPRE -> output is a MatHYPRE */
198: MatMatMult(C,D,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CD);
199: MatMatMult(C,D,MAT_REUSE_MATRIX,PETSC_DEFAULT,&CD);
200: MatAXPY(CD,-1.,pAB,DIFFERENT_NONZERO_PATTERN);
201: MatHasOperation(CD,MATOP_NORM,&flg);
202: if (!flg) { /* TODO add MatNorm_HYPRE */
203: MatConvert(CD,MATAIJ,MAT_INPLACE_MATRIX,&CD);
204: }
205: MatNorm(CD,NORM_INFINITY,&err);
206: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMatMult %g %g",err,norm);
207: MatDestroy(&C);
208: MatDestroy(&D);
209: MatDestroy(&CD);
210: MatDestroy(&pAB);
212: /* When configured with HYPRE, MatMatMatMult is available for the triplet transpose(aij)-aij-aij */
213: MatCreateTranspose(A,&C);
214: MatMatMatMult(C,A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CAB);
215: MatDestroy(&C);
216: MatTranspose(A,MAT_INITIAL_MATRIX,&C);
217: MatMatMult(C,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
218: MatDestroy(&C);
219: MatMatMult(D,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
220: MatNorm(C,NORM_INFINITY,&norm);
221: MatAXPY(C,-1.,CAB,DIFFERENT_NONZERO_PATTERN);
222: MatNorm(C,NORM_INFINITY,&err);
223: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMatMatMult %g %g",err,norm);
224: MatDestroy(&C);
225: MatDestroy(&D);
226: MatDestroy(&CAB);
227: MatDestroy(&B);
229: /* Check MatView */
230: MatViewFromOptions(A,NULL,"-view_A");
231: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
232: MatViewFromOptions(B,NULL,"-view_B");
234: /* Check MatDuplicate/MatCopy */
235: for (j=0;j<3;j++) {
236: MatDuplicateOption dop;
238: dop = MAT_COPY_VALUES;
239: if (j==1) dop = MAT_DO_NOT_COPY_VALUES;
240: if (j==2) dop = MAT_SHARE_NONZERO_PATTERN;
242: for (i=0;i<3;i++) {
243: MatStructure str;
245: PetscPrintf(PETSC_COMM_WORLD,"Dup/Copy tests: %D %D\n",j,i);
247: str = DIFFERENT_NONZERO_PATTERN;
248: if (i==1) str = SAME_NONZERO_PATTERN;
249: if (i==2) str = SUBSET_NONZERO_PATTERN;
251: MatDuplicate(A,dop,&C);
252: MatDuplicate(B,dop,&D);
253: if (dop != MAT_COPY_VALUES) {
254: MatCopy(A,C,str);
255: MatCopy(B,D,str);
256: }
257: /* AXPY with AIJ and HYPRE */
258: MatAXPY(C,-1.0,D,str);
259: MatNorm(C,NORM_INFINITY,&err);
260: if (err > PETSC_SMALL) {
261: MatViewFromOptions(A,NULL,"-view_duplicate_diff");
262: MatViewFromOptions(B,NULL,"-view_duplicate_diff");
263: MatViewFromOptions(C,NULL,"-view_duplicate_diff");
264: SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 1 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
265: }
266: /* AXPY with HYPRE and HYPRE */
267: MatAXPY(D,-1.0,B,str);
268: if (err > PETSC_SMALL) {
269: MatViewFromOptions(A,NULL,"-view_duplicate_diff");
270: MatViewFromOptions(B,NULL,"-view_duplicate_diff");
271: MatViewFromOptions(D,NULL,"-view_duplicate_diff");
272: SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 2 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
273: }
274: /* Copy from HYPRE to AIJ */
275: MatCopy(B,C,str);
276: /* Copy from AIJ to HYPRE */
277: MatCopy(A,D,str);
278: /* AXPY with HYPRE and AIJ */
279: MatAXPY(D,-1.0,C,str);
280: MatHasOperation(D,MATOP_NORM,&flg);
281: if (!flg) { /* TODO add MatNorm_HYPRE */
282: MatConvert(D,MATAIJ,MAT_INPLACE_MATRIX,&D);
283: }
284: MatNorm(D,NORM_INFINITY,&err);
285: if (err > PETSC_SMALL) {
286: MatViewFromOptions(A,NULL,"-view_duplicate_diff");
287: MatViewFromOptions(C,NULL,"-view_duplicate_diff");
288: MatViewFromOptions(D,NULL,"-view_duplicate_diff");
289: SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 3 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
290: }
291: MatDestroy(&C);
292: MatDestroy(&D);
293: }
294: }
295: MatDestroy(&B);
297: MatHasCongruentLayouts(A,&flg);
298: if (flg) {
299: Vec y,y2;
301: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
302: MatCreateVecs(A,NULL,&y);
303: MatCreateVecs(B,NULL,&y2);
304: MatGetDiagonal(A,y);
305: MatGetDiagonal(B,y2);
306: VecAXPY(y2,-1.0,y);
307: VecNorm(y2,NORM_INFINITY,&err);
308: if (err > PETSC_SMALL) {
309: VecViewFromOptions(y,NULL,"-view_diagonal_diff");
310: VecViewFromOptions(y2,NULL,"-view_diagonal_diff");
311: SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatGetDiagonal %g",err);
312: }
313: MatDestroy(&B);
314: VecDestroy(&y);
315: VecDestroy(&y2);
316: }
318: MatDestroy(&A);
320: PetscFinalize();
321: return ierr;
322: }
325: /*TEST
327: build:
328: requires: hypre
330: test:
331: suffix: 1
332: args: -N 11 -M 11
333: output_file: output/ex115_1.out
335: test:
336: suffix: 2
337: nsize: 3
338: args: -N 13 -M 13 -matmatmult_via hypre
339: output_file: output/ex115_1.out
341: test:
342: suffix: 3
343: nsize: 4
344: args: -M 13 -N 7 -matmatmult_via hypre
345: output_file: output/ex115_1.out
347: test:
348: suffix: 4
349: nsize: 2
350: args: -M 12 -N 19
351: output_file: output/ex115_1.out
353: test:
354: suffix: 5
355: nsize: 3
356: args: -M 13 -N 13 -matptap_via hypre -matptap_hypre_outtype hypre
357: output_file: output/ex115_1.out
359: test:
360: suffix: 6
361: nsize: 3
362: args: -M 12 -N 19 -test_offproc
363: output_file: output/ex115_1.out
365: test:
366: suffix: 7
367: nsize: 3
368: args: -M 19 -N 12 -test_offproc -view_B ::ascii_info_detail
369: output_file: output/ex115_7.out
371: test:
372: suffix: 8
373: nsize: 3
374: args: -M 1 -N 12 -test_offproc
375: output_file: output/ex115_1.out
377: test:
378: suffix: 9
379: nsize: 3
380: args: -M 1 -N 2 -test_offproc
381: output_file: output/ex115_1.out
383: TEST*/