Actual source code: ex115.c
petsc-3.10.5 2019-03-28
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: Vec x,x2,y,y2;
12: PetscReal err;
13: PetscInt i,j,N = 6, M = 6;
14: PetscErrorCode ierr;
15: PetscBool fileflg;
16: PetscReal norm;
17: char file[256];
19: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20: PetscOptionsGetString(NULL,NULL,"-f",file,256,&fileflg);
21: MatCreate(PETSC_COMM_WORLD,&A);
22: if (!fileflg) { /* Create a matrix and test MatSetValues */
23: PetscMPIInt NP;
24: PetscBool test_rectangular = PETSC_FALSE;
26: MPI_Comm_size(PETSC_COMM_WORLD,&NP);
27: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
28: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
29: if (N < 6) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Matrix has to have more than 6 columns");
30: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
31: MatSetType(A,MATAIJ);
32: MatSeqAIJSetPreallocation(A,9,NULL);
33: MatMPIAIJSetPreallocation(A,9,NULL,9,NULL);
34: MatCreate(PETSC_COMM_WORLD,&B);
35: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);
36: MatSetType(B,MATHYPRE);
37: if (M == N) {
38: MatHYPRESetPreallocation(B,9,NULL,9,NULL);
39: } else {
40: MatHYPRESetPreallocation(B,6,NULL,6,NULL);
41: }
42: if (M == N) {
43: for (i=0; i<M; i++) {
44: PetscInt cols[] = {0,1,2,3,4,5};
45: PetscScalar vals[] = {0,1./NP,2./NP,3./NP,4./NP,5./NP};
46: for (j=i-2; j<i+1; j++) {
47: if (j >= N) {
48: MatSetValue(A,i,N-1,(1.*j*N+i)/(3.*N*NP),ADD_VALUES);
49: MatSetValue(B,i,N-1,(1.*j*N+i)/(3.*N*NP),ADD_VALUES);
50: } else if (i > j) {
51: MatSetValue(A,i,j,(1.*j*N+i)/(2.*N*NP),ADD_VALUES);
52: MatSetValue(B,i,j,(1.*j*N+i)/(2.*N*NP),ADD_VALUES);
53: } else {
54: MatSetValue(A,i,j,-1.-(1.*j*N+i)/(4.*N*NP),ADD_VALUES);
55: MatSetValue(B,i,j,-1.-(1.*j*N+i)/(4.*N*NP),ADD_VALUES);
56: }
57: }
58: MatSetValues(A,1,&i,6,cols,vals,ADD_VALUES);
59: MatSetValues(B,1,&i,6,cols,vals,ADD_VALUES);
60: }
61: } else { /* HYPRE_IJMatrix does not support INSERT_VALUES with off-proc entries */
62: PetscInt rows[2];
63: MatGetOwnershipRange(A,&rows[0],&rows[1]);
64: for (i=rows[0];i<rows[1];i++) {
65: PetscInt cols[] = {0,1,2,3,4,5};
66: PetscScalar vals[] = {-1,1,-2,2,-3,3};
68: MatSetValues(A,1,&i,6,cols,vals,INSERT_VALUES);
69: MatSetValues(B,1,&i,6,cols,vals,INSERT_VALUES);
70: }
71: }
72: /* MAT_FLUSH_ASSEMBLY currently not supported */
73: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
74: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
75: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
78: PetscOptionsGetBool(NULL,NULL,"-test_rectangular",&test_rectangular,NULL);
79: if (M != N && !test_rectangular) {
80: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
81: MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
82: } else {
83: /* MatAXPY_Basic further exercises MatSetValues_HYPRE */
84: /* however there are still some issues
85: with rectangular matrices.
86: try, e.g., mpiexec -n 3 ./ex115 -M 7 -N 6 -test_rectangular */
87: MatAXPY(B,-1.,A,DIFFERENT_NONZERO_PATTERN);
88: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
89: }
90: MatNorm(C,NORM_INFINITY,&err);
91: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatSetValues %g",err);
92: MatDestroy(&B);
93: MatDestroy(&C);
94: } else {
95: PetscViewer viewer;
96: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&viewer);
97: MatSetFromOptions(A);
98: MatLoad(A,viewer);
99: PetscViewerDestroy(&viewer);
100: MatGetSize(A,&M,&N);
101: }
102: /* check conversion routines */
103: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&B);
104: MatConvert(A,MATHYPRE,MAT_REUSE_MATRIX,&B);
105: MatConvert(B,MATIS,MAT_INITIAL_MATRIX,&D);
106: MatConvert(B,MATIS,MAT_REUSE_MATRIX,&D);
107: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&C);
108: MatConvert(B,MATAIJ,MAT_REUSE_MATRIX,&C);
109: MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
110: MatNorm(C,NORM_INFINITY,&err);
111: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat AIJ %g",err);
112: MatDestroy(&C);
113: MatConvert(D,MATAIJ,MAT_INITIAL_MATRIX,&C);
114: MatAXPY(C,-1.,A,SAME_NONZERO_PATTERN);
115: MatNorm(C,NORM_INFINITY,&err);
116: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error Mat IS %g",err);
117: MatDestroy(&C);
118: MatDestroy(&D);
120: /* check MatCreateFromParCSR */
121: MatHYPREGetParCSR(B,&parcsr);
122: MatCreateFromParCSR(parcsr,MATAIJ,PETSC_COPY_VALUES,&D);
123: MatDestroy(&D);
124: MatCreateFromParCSR(parcsr,MATHYPRE,PETSC_USE_POINTER,&C);
126: /* check matmult */
127: MatCreateVecs(A,&x,&y);
128: MatCreateVecs(A,&x2,&y2);
129: VecSet(x,1.);
130: MatMult(A,x,y);
131: MatMult(B,x,y2);
132: VecAXPY(y,-1.,y2);
133: VecNorm(y,NORM_2,&err);
134: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMult B %g",err);
135: MatMult(C,x,y);
136: VecAXPY(y,-1.,y2);
137: VecNorm(y,NORM_2,&err);
138: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMult C %g",err);
139: VecSet(y,1.);
140: MatMultTranspose(A,y,x);
141: MatMultTranspose(B,y,x2);
142: VecAXPY(x,-1.,x2);
143: VecNorm(x,NORM_2,&err);
144: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTranspose C %g",err);
145: MatMultTranspose(C,y,x);
146: VecAXPY(x,-1.,x2);
147: VecNorm(x,NORM_2,&err);
148: if (err > PETSC_SMALL) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMultTranspose C %g",err);
150: /* check PtAP */
151: if (M == N) {
152: Mat pP,hP;
154: /* PETSc MatPtAP -> output is a MatAIJ
155: It uses HYPRE functions when -matptap_via hypre is specified at command line */
156: MatPtAP(A,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pP);
157: MatPtAP(A,A,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pP);
158: MatNorm(pP,NORM_INFINITY,&norm);
160: /* MatPtAP_HYPRE_HYPRE -> output is a MatHYPRE */
161: MatPtAP(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
162: MatPtAP(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
163: MatConvert(hP,MATAIJ,MAT_INITIAL_MATRIX,&D);
164: MatAXPY(D,-1.,pP,DIFFERENT_NONZERO_PATTERN);
165: MatNorm(D,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);
168: MatDestroy(&D);
170: /* MatPtAP_AIJ_HYPRE -> output can be decided at runtime with -matptap_hypre_outtype */
171: MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&hP);
172: MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&hP);
173: MatConvert(hP,MATAIJ,MAT_INITIAL_MATRIX,&D);
174: MatAXPY(D,-1.,pP,DIFFERENT_NONZERO_PATTERN);
175: MatNorm(D,NORM_INFINITY,&err);
176: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatPtAP mixed %g %g",err,norm);
177: MatDestroy(&D);
178: MatDestroy(&hP);
180: MatDestroy(&pP);
181: }
183: /* check MatMatMult */
184: MatDestroy(&C);
185: MatDestroy(&B);
186: if (fileflg) {
187: PetscOptionsGetString(NULL,NULL,"-fB",file,256,&fileflg);
188: if (fileflg) {
189: PetscViewer viewer;
190: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&viewer);
191: MatSetFromOptions(B);
192: MatLoad(B,viewer);
193: PetscViewerDestroy(&viewer);
194: }
195: }
196: if (!B) {
197: MatTranspose(A,MAT_INITIAL_MATRIX,&B);
198: }
199: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&C);
200: MatConvert(B,MATHYPRE,MAT_INITIAL_MATRIX,&D);
202: /* PETSc MatMatMult -> output is a MatAIJ
203: It uses HYPRE functions when -matmatmult_via hypre is specified at command line */
204: MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pAB);
205: MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pAB);
206: MatNorm(pAB,NORM_INFINITY,&norm);
208: /* MatMatMult_HYPRE_HYPRE -> output is a MatHYPRE */
209: MatMatMult(C,D,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CD);
210: MatMatMult(C,D,MAT_REUSE_MATRIX,PETSC_DEFAULT,&CD);
211: MatDestroy(&D);
212: MatConvert(CD,MATAIJ,MAT_INITIAL_MATRIX,&D);
213: MatAXPY(D,-1.,pAB,DIFFERENT_NONZERO_PATTERN);
214: MatNorm(D,NORM_INFINITY,&err);
215: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMatMult %g %g",err,norm);
216: MatDestroy(&C);
217: MatDestroy(&D);
218: MatDestroy(&CD);
219: MatDestroy(&pAB);
221: /* When configured with HYPRE, MatMatMatMult is available for the triplet transpose(aij)-aij-aij */
222: MatCreateTranspose(A,&C);
223: MatMatMatMult(C,A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CAB);
224: MatDestroy(&C);
225: MatTranspose(A,MAT_INITIAL_MATRIX,&C);
226: MatMatMult(C,A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
227: MatDestroy(&C);
228: MatMatMult(D,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
229: MatNorm(C,NORM_INFINITY,&norm);
230: MatAXPY(C,-1.,CAB,DIFFERENT_NONZERO_PATTERN);
231: MatNorm(C,NORM_INFINITY,&err);
232: if (err/norm > PETSC_SMALL) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Error MatMatMatMult %g %g",err,norm);
233: MatDestroy(&C);
234: MatDestroy(&D);
235: MatDestroy(&CAB);
237: VecDestroy(&x);
238: VecDestroy(&x2);
239: VecDestroy(&y);
240: VecDestroy(&y2);
241: MatDestroy(&A);
242: MatDestroy(&B);
244: PetscFinalize();
245: return ierr;
246: }
249: /*TEST
251: build:
252: requires: hypre
254: test:
255: suffix: 1
256: requires: hypre
257: args: -N 11 -M 11
259: test:
260: suffix: 2
261: nsize: 3
262: requires: hypre
263: args: -N 13 -M 13 -matmatmult_via hypre
265: test:
266: suffix: 3
267: nsize: 4
268: requires: hypre
269: args: -M 13 -N 7 -matmatmult_via hypre
271: test:
272: suffix: 4
273: nsize: 2
274: requires: hypre
275: args: -M 12 -N 19
277: test:
278: suffix: 5
279: nsize: 3
280: requires: hypre
281: args: -M 13 -N 13 -matptap_via hypre -matptap_hypre_outtype hypre
283: TEST*/