Actual source code: ex115.c
petsc-3.14.6 2021-03-30
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: 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);
169: MatPtAPMultEqual(A,A,pP,10,&flg);
170: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP_MatAIJ");
172: /* MatPtAP_HYPRE_HYPRE -> output is a MatHYPRE */
173: MatPtAP(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
174: MatPtAP(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
175: MatPtAPMultEqual(C,B,hP,10,&flg);
176: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP_HYPRE_HYPRE");
178: /* Test MatAXPY_Basic() */
179: MatAXPY(hP,-1.,pP,DIFFERENT_NONZERO_PATTERN);
180: MatHasOperation(hP,MATOP_NORM,&flg);
181: if (!flg) { /* TODO add MatNorm_HYPRE */
182: MatConvert(hP,MATAIJ,MAT_INPLACE_MATRIX,&hP);
183: }
184: MatNorm(hP,NORM_INFINITY,&err);
185: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)hP),PETSC_ERR_PLIB,"Error MatPtAP %g %g",err,norm);
186: MatDestroy(&pP);
187: MatDestroy(&hP);
189: /* MatPtAP_AIJ_HYPRE -> output can be decided at runtime with -matptap_hypre_outtype */
190: MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
191: MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
192: MatPtAPMultEqual(A,B,hP,10,&flg);
193: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP_AIJ_HYPRE");
194: MatDestroy(&hP);
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);
210: MatMatMultEqual(A,B,pAB,10,&flg);
211: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatMult_AIJ_AIJ");
213: /* MatMatMult_HYPRE_HYPRE -> output is a MatHYPRE */
214: MatMatMult(C,D,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CD);
215: MatMatMult(C,D,MAT_REUSE_MATRIX,PETSC_DEFAULT,&CD);
216: MatMatMultEqual(C,D,CD,10,&flg);
217: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatMult_HYPRE_HYPRE");
219: /* Test MatAXPY_Basic() */
220: MatAXPY(CD,-1.,pAB,DIFFERENT_NONZERO_PATTERN);
222: MatHasOperation(CD,MATOP_NORM,&flg);
223: if (!flg) { /* TODO add MatNorm_HYPRE */
224: MatConvert(CD,MATAIJ,MAT_INPLACE_MATRIX,&CD);
225: }
226: MatNorm(CD,NORM_INFINITY,&err);
227: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)CD),PETSC_ERR_PLIB,"Error MatMatMult %g %g",err,norm);
229: MatDestroy(&C);
230: MatDestroy(&D);
231: MatDestroy(&pAB);
232: MatDestroy(&CD);
234: /* When configured with HYPRE, MatMatMatMult is available for the triplet transpose(aij)-aij-aij */
235: MatCreateTranspose(A,&C);
236: MatMatMatMult(C,A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CAB);
237: MatDestroy(&C);
238: MatTranspose(A,MAT_INITIAL_MATRIX,&C);
239: MatMatMult(C,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
240: MatDestroy(&C);
241: MatMatMult(D,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
242: MatNorm(C,NORM_INFINITY,&norm);
243: MatAXPY(C,-1.,CAB,DIFFERENT_NONZERO_PATTERN);
244: MatNorm(C,NORM_INFINITY,&err);
245: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMatMatMult %g %g",err,norm);
246: MatDestroy(&C);
247: MatDestroy(&D);
248: MatDestroy(&CAB);
249: MatDestroy(&B);
250: }
252: /* Check MatView */
253: MatViewFromOptions(A,NULL,"-view_A");
254: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
255: MatViewFromOptions(B,NULL,"-view_B");
257: /* Check MatDuplicate/MatCopy */
258: for (j=0;j<3;j++) {
259: MatDuplicateOption dop;
261: dop = MAT_COPY_VALUES;
262: if (j==1) dop = MAT_DO_NOT_COPY_VALUES;
263: if (j==2) dop = MAT_SHARE_NONZERO_PATTERN;
265: for (i=0;i<3;i++) {
266: MatStructure str;
268: PetscPrintf(PETSC_COMM_WORLD,"Dup/Copy tests: %D %D\n",j,i);
270: str = DIFFERENT_NONZERO_PATTERN;
271: if (i==1) str = SAME_NONZERO_PATTERN;
272: if (i==2) str = SUBSET_NONZERO_PATTERN;
274: MatDuplicate(A,dop,&C);
275: MatDuplicate(B,dop,&D);
276: if (dop != MAT_COPY_VALUES) {
277: MatCopy(A,C,str);
278: MatCopy(B,D,str);
279: }
280: /* AXPY with AIJ and HYPRE */
281: MatAXPY(C,-1.0,D,str);
282: MatNorm(C,NORM_INFINITY,&err);
283: if (err > PETSC_SMALL) {
284: MatViewFromOptions(A,NULL,"-view_duplicate_diff");
285: MatViewFromOptions(B,NULL,"-view_duplicate_diff");
286: MatViewFromOptions(C,NULL,"-view_duplicate_diff");
287: SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 1 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
288: }
289: /* AXPY with HYPRE and HYPRE */
290: MatAXPY(D,-1.0,B,str);
291: if (err > PETSC_SMALL) {
292: MatViewFromOptions(A,NULL,"-view_duplicate_diff");
293: MatViewFromOptions(B,NULL,"-view_duplicate_diff");
294: MatViewFromOptions(D,NULL,"-view_duplicate_diff");
295: SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 2 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
296: }
297: /* Copy from HYPRE to AIJ */
298: MatCopy(B,C,str);
299: /* Copy from AIJ to HYPRE */
300: MatCopy(A,D,str);
301: /* AXPY with HYPRE and AIJ */
302: MatAXPY(D,-1.0,C,str);
303: MatHasOperation(D,MATOP_NORM,&flg);
304: if (!flg) { /* TODO add MatNorm_HYPRE */
305: MatConvert(D,MATAIJ,MAT_INPLACE_MATRIX,&D);
306: }
307: MatNorm(D,NORM_INFINITY,&err);
308: if (err > PETSC_SMALL) {
309: MatViewFromOptions(A,NULL,"-view_duplicate_diff");
310: MatViewFromOptions(C,NULL,"-view_duplicate_diff");
311: MatViewFromOptions(D,NULL,"-view_duplicate_diff");
312: SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error test 3 MatDuplicate/MatCopy %g (%D,%D)",err,j,i);
313: }
314: MatDestroy(&C);
315: MatDestroy(&D);
316: }
317: }
318: MatDestroy(&B);
320: MatHasCongruentLayouts(A,&flg);
321: if (flg) {
322: Vec y,y2;
324: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
325: MatCreateVecs(A,NULL,&y);
326: MatCreateVecs(B,NULL,&y2);
327: MatGetDiagonal(A,y);
328: MatGetDiagonal(B,y2);
329: VecAXPY(y2,-1.0,y);
330: VecNorm(y2,NORM_INFINITY,&err);
331: if (err > PETSC_SMALL) {
332: VecViewFromOptions(y,NULL,"-view_diagonal_diff");
333: VecViewFromOptions(y2,NULL,"-view_diagonal_diff");
334: SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatGetDiagonal %g",err);
335: }
336: MatDestroy(&B);
337: VecDestroy(&y);
338: VecDestroy(&y2);
339: }
341: MatDestroy(&A);
343: PetscFinalize();
344: return ierr;
345: }
348: /*TEST
350: build:
351: requires: hypre
353: test:
354: suffix: 1
355: args: -N 11 -M 11
356: output_file: output/ex115_1.out
358: test:
359: suffix: 2
360: nsize: 3
361: args: -N 13 -M 13 -matmatmult_via hypre
362: output_file: output/ex115_1.out
364: test:
365: suffix: 3
366: nsize: 4
367: args: -M 13 -N 7 -matmatmult_via hypre
368: output_file: output/ex115_1.out
370: test:
371: suffix: 4
372: nsize: 2
373: args: -M 12 -N 19
374: output_file: output/ex115_1.out
376: test:
377: suffix: 5
378: nsize: 3
379: args: -M 13 -N 13 -matptap_via hypre -matptap_hypre_outtype hypre
380: output_file: output/ex115_1.out
382: test:
383: suffix: 6
384: nsize: 3
385: args: -M 12 -N 19 -test_offproc
386: output_file: output/ex115_1.out
388: test:
389: suffix: 7
390: nsize: 3
391: args: -M 19 -N 12 -test_offproc -view_B ::ascii_info_detail
392: output_file: output/ex115_7.out
394: test:
395: suffix: 8
396: nsize: 3
397: args: -M 1 -N 12 -test_offproc
398: output_file: output/ex115_1.out
400: test:
401: suffix: 9
402: nsize: 3
403: args: -M 1 -N 2 -test_offproc
404: output_file: output/ex115_1.out
406: TEST*/