Actual source code: ex96.c
petsc-3.13.6 2020-09-29
2: static char help[] ="Tests sequential and parallel DMCreateMatrix(), MatMatMult() and MatPtAP()\n\
3: -Mx <xg>, where <xg> = number of coarse grid points in the x-direction\n\
4: -My <yg>, where <yg> = number of coarse grid points in the y-direction\n\
5: -Mz <zg>, where <zg> = number of coarse grid points in the z-direction\n\
6: -Npx <npx>, where <npx> = number of processors in the x-direction\n\
7: -Npy <npy>, where <npy> = number of processors in the y-direction\n\
8: -Npz <npz>, where <npz> = number of processors in the z-direction\n\n";
10: /*
11: This test is modified from ~src/ksp/tests/ex19.c.
12: Example of usage: mpiexec -n 3 ./ex96 -Mx 10 -My 10 -Mz 10
13: */
15: #include <petscdm.h>
16: #include <petscdmda.h>
18: /* User-defined Section 1.5 Writing Application Codes with PETSc contexts */
19: typedef struct {
20: PetscInt mx,my,mz; /* number grid points in x, y and z direction */
21: Vec localX,localF; /* local vectors with ghost region */
22: DM da;
23: Vec x,b,r; /* global vectors */
24: Mat J; /* Jacobian on grid */
25: } GridCtx;
26: typedef struct {
27: GridCtx fine;
28: GridCtx coarse;
29: PetscInt ratio;
30: Mat Ii; /* interpolation from coarse to fine */
31: } AppCtx;
33: #define COARSE_LEVEL 0
34: #define FINE_LEVEL 1
36: /*
37: Mm_ratio - ration of grid lines between fine and coarse grids.
38: */
39: int main(int argc,char **argv)
40: {
42: AppCtx user;
43: PetscInt Npx=PETSC_DECIDE,Npy=PETSC_DECIDE,Npz=PETSC_DECIDE;
44: PetscMPIInt size,rank;
45: PetscInt m,n,M,N,i,nrows;
46: PetscScalar one = 1.0;
47: PetscReal fill=2.0;
48: Mat A,A_tmp,P,C,C1,C2;
49: PetscScalar *array,none = -1.0,alpha;
50: Vec x,v1,v2,v3,v4;
51: PetscReal norm,norm_tmp,norm_tmp1,tol=100.*PETSC_MACHINE_EPSILON;
52: PetscRandom rdm;
53: PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_3D=PETSC_TRUE,flg;
54: const PetscInt *ia,*ja;
56: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
57: PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);
59: user.ratio = 2;
60: user.coarse.mx = 20; user.coarse.my = 20; user.coarse.mz = 20;
62: PetscOptionsGetInt(NULL,NULL,"-Mx",&user.coarse.mx,NULL);
63: PetscOptionsGetInt(NULL,NULL,"-My",&user.coarse.my,NULL);
64: PetscOptionsGetInt(NULL,NULL,"-Mz",&user.coarse.mz,NULL);
65: PetscOptionsGetInt(NULL,NULL,"-ratio",&user.ratio,NULL);
67: if (user.coarse.mz) Test_3D = PETSC_TRUE;
69: MPI_Comm_size(PETSC_COMM_WORLD,&size);
70: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
71: PetscOptionsGetInt(NULL,NULL,"-Npx",&Npx,NULL);
72: PetscOptionsGetInt(NULL,NULL,"-Npy",&Npy,NULL);
73: PetscOptionsGetInt(NULL,NULL,"-Npz",&Npz,NULL);
75: /* Set up distributed array for fine grid */
76: if (!Test_3D) {
77: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,user.coarse.my,Npx,Npy,1,1,NULL,NULL,&user.coarse.da);
78: } else {
79: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,user.coarse.my,user.coarse.mz,Npx,Npy,Npz,1,1,NULL,NULL,NULL,&user.coarse.da);
80: }
81: DMSetFromOptions(user.coarse.da);
82: DMSetUp(user.coarse.da);
84: /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
85: DMRefine(user.coarse.da,PetscObjectComm((PetscObject)user.coarse.da),&user.fine.da);
87: /* Test DMCreateMatrix() */
88: /*------------------------------------------------------------*/
89: DMSetMatType(user.fine.da,MATAIJ);
90: DMCreateMatrix(user.fine.da,&A);
91: DMSetMatType(user.fine.da,MATBAIJ);
92: DMCreateMatrix(user.fine.da,&C);
94: MatConvert(C,MATAIJ,MAT_INITIAL_MATRIX,&A_tmp); /* not work for mpisbaij matrix! */
95: MatEqual(A,A_tmp,&flg);
96: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"A != C");
97: MatDestroy(&C);
98: MatDestroy(&A_tmp);
100: /*------------------------------------------------------------*/
102: MatGetLocalSize(A,&m,&n);
103: MatGetSize(A,&M,&N);
104: /* if (!rank) printf("A %d, %d\n",M,N); */
106: /* set val=one to A */
107: if (size == 1) {
108: MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
109: if (flg) {
110: MatSeqAIJGetArray(A,&array);
111: for (i=0; i<ia[nrows]; i++) array[i] = one;
112: MatSeqAIJRestoreArray(A,&array);
113: }
114: MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
115: } else {
116: Mat AA,AB;
117: MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);
118: MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
119: if (flg) {
120: MatSeqAIJGetArray(AA,&array);
121: for (i=0; i<ia[nrows]; i++) array[i] = one;
122: MatSeqAIJRestoreArray(AA,&array);
123: }
124: MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
125: MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
126: if (flg) {
127: MatSeqAIJGetArray(AB,&array);
128: for (i=0; i<ia[nrows]; i++) array[i] = one;
129: MatSeqAIJRestoreArray(AB,&array);
130: }
131: MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
132: }
133: /* MatView(A, PETSC_VIEWER_STDOUT_WORLD); */
135: /* Create interpolation between the fine and coarse grids */
136: DMCreateInterpolation(user.coarse.da,user.fine.da,&P,NULL);
137: MatGetLocalSize(P,&m,&n);
138: MatGetSize(P,&M,&N);
139: /* if (!rank) printf("P %d, %d\n",M,N); */
141: /* Create vectors v1 and v2 that are compatible with A */
142: VecCreate(PETSC_COMM_WORLD,&v1);
143: MatGetLocalSize(A,&m,NULL);
144: VecSetSizes(v1,m,PETSC_DECIDE);
145: VecSetFromOptions(v1);
146: VecDuplicate(v1,&v2);
147: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
148: PetscRandomSetFromOptions(rdm);
150: /* Test MatMatMult(): C = A*P */
151: /*----------------------------*/
152: if (Test_MatMatMult) {
153: MatDuplicate(A,MAT_COPY_VALUES,&A_tmp);
154: MatMatMult(A_tmp,P,MAT_INITIAL_MATRIX,fill,&C);
156: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
157: alpha=1.0;
158: for (i=0; i<2; i++) {
159: alpha -= 0.1;
160: MatScale(A_tmp,alpha);
161: MatMatMult(A_tmp,P,MAT_REUSE_MATRIX,fill,&C);
162: }
163: /* Free intermediate data structures created for reuse of C=Pt*A*P */
164: MatFreeIntermediateDataStructures(C);
166: /* Test MatDuplicate() */
167: /*----------------------------*/
168: MatDuplicate(C,MAT_COPY_VALUES,&C1);
169: MatDuplicate(C1,MAT_COPY_VALUES,&C2);
170: MatDestroy(&C1);
171: MatDestroy(&C2);
173: /* Create vector x that is compatible with P */
174: VecCreate(PETSC_COMM_WORLD,&x);
175: MatGetLocalSize(P,NULL,&n);
176: VecSetSizes(x,n,PETSC_DECIDE);
177: VecSetFromOptions(x);
179: norm = 0.0;
180: for (i=0; i<10; i++) {
181: VecSetRandom(x,rdm);
182: MatMult(P,x,v1);
183: MatMult(A_tmp,v1,v2); /* v2 = A*P*x */
184: MatMult(C,x,v1); /* v1 = C*x */
185: VecAXPY(v1,none,v2);
186: VecNorm(v1,NORM_1,&norm_tmp);
187: VecNorm(v2,NORM_1,&norm_tmp1);
188: norm_tmp /= norm_tmp1;
189: if (norm_tmp > norm) norm = norm_tmp;
190: }
191: if (norm >= tol && !rank) {
192: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|/|v2|: %g\n",(double)norm);
193: }
195: VecDestroy(&x);
196: MatDestroy(&C);
197: MatDestroy(&A_tmp);
198: }
200: /* Test P^T * A * P - MatPtAP() */
201: /*------------------------------*/
202: if (Test_MatPtAP) {
203: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
204: MatGetLocalSize(C,&m,&n);
206: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
207: alpha=1.0;
208: for (i=0; i<1; i++) {
209: alpha -= 0.1;
210: MatScale(A,alpha);
211: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
212: }
214: /* Free intermediate data structures created for reuse of C=Pt*A*P */
215: MatFreeIntermediateDataStructures(C);
217: /* Test MatDuplicate() */
218: /*----------------------------*/
219: MatDuplicate(C,MAT_COPY_VALUES,&C1);
220: MatDuplicate(C1,MAT_COPY_VALUES,&C2);
221: MatDestroy(&C1);
222: MatDestroy(&C2);
224: /* Create vector x that is compatible with P */
225: VecCreate(PETSC_COMM_WORLD,&x);
226: MatGetLocalSize(P,&m,&n);
227: VecSetSizes(x,n,PETSC_DECIDE);
228: VecSetFromOptions(x);
230: VecCreate(PETSC_COMM_WORLD,&v3);
231: VecSetSizes(v3,n,PETSC_DECIDE);
232: VecSetFromOptions(v3);
233: VecDuplicate(v3,&v4);
235: norm = 0.0;
236: for (i=0; i<10; i++) {
237: VecSetRandom(x,rdm);
238: MatMult(P,x,v1);
239: MatMult(A,v1,v2); /* v2 = A*P*x */
241: MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
242: MatMult(C,x,v4); /* v3 = C*x */
243: VecAXPY(v4,none,v3);
244: VecNorm(v4,NORM_1,&norm_tmp);
245: VecNorm(v3,NORM_1,&norm_tmp1);
247: norm_tmp /= norm_tmp1;
248: if (norm_tmp > norm) norm = norm_tmp;
249: }
250: if (norm >= tol && !rank) {
251: PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v3 - v4|/|v3|: %g\n",(double)norm);
252: }
253: MatDestroy(&C);
254: VecDestroy(&v3);
255: VecDestroy(&v4);
256: VecDestroy(&x);
257: }
259: /* Clean up */
260: MatDestroy(&A);
261: PetscRandomDestroy(&rdm);
262: VecDestroy(&v1);
263: VecDestroy(&v2);
264: DMDestroy(&user.fine.da);
265: DMDestroy(&user.coarse.da);
266: MatDestroy(&P);
267: PetscFinalize();
268: return ierr;
269: }
271: /*TEST
273: test:
274: args: -Mx 10 -My 5 -Mz 10
275: output_file: output/ex96_1.out
277: test:
278: suffix: nonscalable
279: nsize: 3
280: args: -Mx 10 -My 5 -Mz 10
281: output_file: output/ex96_1.out
283: test:
284: suffix: scalable
285: nsize: 3
286: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable
287: output_file: output/ex96_1.out
289: test:
290: suffix: seq_scalable
291: nsize: 3
292: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via scalable -inner_offdiag_matproduct_ab_via scalable
293: output_file: output/ex96_1.out
295: test:
296: suffix: seq_sorted
297: nsize: 3
298: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via sorted -inner_offdiag_matproduct_ab_via sorted
299: output_file: output/ex96_1.out
301: test:
302: suffix: seq_scalable_fast
303: nsize: 3
304: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via scalable_fast -inner_offdiag_matproduct_ab_via scalable_fast
305: output_file: output/ex96_1.out
307: test:
308: suffix: seq_heap
309: nsize: 3
310: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via heap -inner_offdiag_matproduct_ab_via heap
311: output_file: output/ex96_1.out
313: test:
314: suffix: seq_btheap
315: nsize: 3
316: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via btheap -inner_offdiag_matproduct_ab_via btheap
317: output_file: output/ex96_1.out
319: test:
320: suffix: seq_llcondensed
321: nsize: 3
322: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via llcondensed -inner_offdiag_matproduct_ab_via llcondensed
323: output_file: output/ex96_1.out
325: test:
326: suffix: seq_rowmerge
327: nsize: 3
328: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via rowmerge -inner_offdiag_matproduct_ab_via rowmerge
329: output_file: output/ex96_1.out
331: test:
332: suffix: allatonce
333: nsize: 3
334: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via allatonce
335: output_file: output/ex96_1.out
337: test:
338: suffix: allatonce_merged
339: nsize: 3
340: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via allatonce_merged
341: output_file: output/ex96_1.out
343: TEST*/