Actual source code: ex96.c
petsc-3.6.4 2016-04-12
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/examples/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>
17: #include <../src/mat/impls/aij/seq/aij.h>
18: #include <../src/mat/impls/aij/mpi/mpiaij.h>
20: /* User-defined application contexts */
21: typedef struct {
22: PetscInt mx,my,mz; /* number grid points in x, y and z direction */
23: Vec localX,localF; /* local vectors with ghost region */
24: DM da;
25: Vec x,b,r; /* global vectors */
26: Mat J; /* Jacobian on grid */
27: } GridCtx;
28: typedef struct {
29: GridCtx fine;
30: GridCtx coarse;
31: PetscInt ratio;
32: Mat Ii; /* interpolation from coarse to fine */
33: } AppCtx;
35: #define COARSE_LEVEL 0
36: #define FINE_LEVEL 1
38: /*
39: Mm_ratio - ration of grid lines between fine and coarse grids.
40: */
43: int main(int argc,char **argv)
44: {
46: AppCtx user;
47: PetscInt Npx=PETSC_DECIDE,Npy=PETSC_DECIDE,Npz=PETSC_DECIDE;
48: PetscMPIInt size,rank;
49: PetscInt m,n,M,N,i,nrows;
50: PetscScalar one = 1.0;
51: PetscReal fill=2.0;
52: Mat A,A_tmp,P,C,C1,C2;
53: PetscScalar *array,none = -1.0,alpha;
54: Vec x,v1,v2,v3,v4;
55: PetscReal norm,norm_tmp,norm_tmp1,tol=1.e-12;
56: PetscRandom rdm;
57: PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_3D=PETSC_FALSE,flg;
59: PetscInitialize(&argc,&argv,NULL,help);
60: PetscOptionsGetReal(NULL,"-tol",&tol,NULL);
62: user.ratio = 2;
63: user.coarse.mx = 2; user.coarse.my = 2; user.coarse.mz = 0;
65: PetscOptionsGetInt(NULL,"-Mx",&user.coarse.mx,NULL);
66: PetscOptionsGetInt(NULL,"-My",&user.coarse.my,NULL);
67: PetscOptionsGetInt(NULL,"-Mz",&user.coarse.mz,NULL);
68: PetscOptionsGetInt(NULL,"-ratio",&user.ratio,NULL);
70: if (user.coarse.mz) Test_3D = PETSC_TRUE;
72: user.fine.mx = user.ratio*(user.coarse.mx-1)+1;
73: user.fine.my = user.ratio*(user.coarse.my-1)+1;
74: user.fine.mz = user.ratio*(user.coarse.mz-1)+1;
76: MPI_Comm_size(PETSC_COMM_WORLD,&size);
77: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
78: PetscOptionsGetInt(NULL,"-Npx",&Npx,NULL);
79: PetscOptionsGetInt(NULL,"-Npy",&Npy,NULL);
80: PetscOptionsGetInt(NULL,"-Npz",&Npz,NULL);
82: /* Set up distributed array for fine grid */
83: if (!Test_3D) {
84: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.fine.mx,
85: user.fine.my,Npx,Npy,1,1,NULL,NULL,&user.fine.da);
86: } else {
87: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,
88: user.fine.mx,user.fine.my,user.fine.mz,Npx,Npy,Npz,
89: 1,1,NULL,NULL,NULL,&user.fine.da);
90: }
92: /* Test DMCreateMatrix() */
93: /*------------------------------------------------------------*/
94: DMSetMatType(user.fine.da,MATAIJ);
95: DMCreateMatrix(user.fine.da,&A);
96: DMSetMatType(user.fine.da,MATBAIJ);
97: DMCreateMatrix(user.fine.da,&C);
99: MatConvert(C,MATAIJ,MAT_INITIAL_MATRIX,&A_tmp); /* not work for mpisbaij matrix! */
100: MatEqual(A,A_tmp,&flg);
101: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"A != C");
102: MatDestroy(&C);
103: MatDestroy(&A_tmp);
105: /*------------------------------------------------------------*/
107: MatGetLocalSize(A,&m,&n);
108: MatGetSize(A,&M,&N);
109: /* set val=one to A */
110: if (size == 1) {
111: const PetscInt *ia,*ja;
112: MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
113: if (flg) {
114: MatSeqAIJGetArray(A,&array);
115: for (i=0; i<ia[nrows]; i++) array[i] = one;
116: MatSeqAIJRestoreArray(A,&array);
117: }
118: MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
119: } else {
120: Mat_MPIAIJ *aij = (Mat_MPIAIJ*)A->data;
121: Mat_SeqAIJ *a = (Mat_SeqAIJ*)(aij->A)->data, *b=(Mat_SeqAIJ*)(aij->B)->data;
122: /* A_part */
123: for (i=0; i<a->i[m]; i++) a->a[i] = one;
124: /* B_part */
125: for (i=0; i<b->i[m]; i++) b->a[i] = one;
127: }
128: /* MatView(A, PETSC_VIEWER_STDOUT_WORLD); */
130: /* Set up distributed array for coarse grid */
131: if (!Test_3D) {
132: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,
133: user.coarse.my,Npx,Npy,1,1,NULL,NULL,&user.coarse.da);
134: } else {
135: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,
136: user.coarse.mx,user.coarse.my,user.coarse.mz,Npx,Npy,Npz,
137: 1,1,NULL,NULL,NULL,&user.coarse.da);
138: }
140: /* Create interpolation between the levels */
141: DMCreateInterpolation(user.coarse.da,user.fine.da,&P,NULL);
143: MatGetLocalSize(P,&m,&n);
144: MatGetSize(P,&M,&N);
146: /* Create vectors v1 and v2 that are compatible with A */
147: VecCreate(PETSC_COMM_WORLD,&v1);
148: MatGetLocalSize(A,&m,NULL);
149: VecSetSizes(v1,m,PETSC_DECIDE);
150: VecSetFromOptions(v1);
151: VecDuplicate(v1,&v2);
152: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
153: PetscRandomSetFromOptions(rdm);
155: /* Test MatMatMult(): C = A*P */
156: /*----------------------------*/
157: if (Test_MatMatMult) {
158: MatDuplicate(A,MAT_COPY_VALUES,&A_tmp);
159: MatMatMult(A_tmp,P,MAT_INITIAL_MATRIX,fill,&C);
161: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
162: alpha=1.0;
163: for (i=0; i<2; i++) {
164: alpha -=0.1;
165: MatScale(A_tmp,alpha);
166: MatMatMult(A_tmp,P,MAT_REUSE_MATRIX,fill,&C);
167: }
169: /* Test MatDuplicate() */
170: /*----------------------------*/
171: MatDuplicate(C,MAT_COPY_VALUES,&C1);
172: MatDuplicate(C1,MAT_COPY_VALUES,&C2);
173: MatDestroy(&C1);
174: MatDestroy(&C2);
176: /* Create vector x that is compatible with P */
177: VecCreate(PETSC_COMM_WORLD,&x);
178: MatGetLocalSize(P,NULL,&n);
179: VecSetSizes(x,n,PETSC_DECIDE);
180: VecSetFromOptions(x);
182: norm = 0.0;
183: for (i=0; i<10; i++) {
184: VecSetRandom(x,rdm);
185: MatMult(P,x,v1);
186: MatMult(A_tmp,v1,v2); /* v2 = A*P*x */
187: MatMult(C,x,v1); /* v1 = C*x */
188: VecAXPY(v1,none,v2);
189: VecNorm(v1,NORM_1,&norm_tmp);
190: VecNorm(v2,NORM_1,&norm_tmp1);
191: norm_tmp /= norm_tmp1;
192: if (norm_tmp > norm) norm = norm_tmp;
193: }
194: if (norm >= tol && !rank) {
195: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|/|v2|: %g\n",(double)norm);
196: }
198: VecDestroy(&x);
199: MatDestroy(&C);
200: MatDestroy(&A_tmp);
201: }
203: /* Test P^T * A * P - MatPtAP() */
204: /*------------------------------*/
205: if (Test_MatPtAP) {
206: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
207: MatGetLocalSize(C,&m,&n);
209: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
210: alpha=1.0;
211: for (i=0; i<1; i++) {
212: alpha -=0.1;
213: MatScale(A,alpha);
214: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
215: }
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 0;
269: }