Actual source code: ex96.c

petsc-3.4.2 2013-07-02
  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 <petscdmda.h>
 16: #include <../src/mat/impls/aij/seq/aij.h>
 17: #include <../src/mat/impls/aij/mpi/mpiaij.h>

 19: /* User-defined application contexts */
 20: typedef struct {
 21:   PetscInt mx,my,mz;            /* number grid points in x, y and z direction */
 22:   Vec      localX,localF;       /* local vectors with ghost region */
 23:   DM       da;
 24:   Vec      x,b,r;               /* global vectors */
 25:   Mat      J;                   /* Jacobian on grid */
 26: } GridCtx;
 27: typedef struct {
 28:   GridCtx  fine;
 29:   GridCtx  coarse;
 30:   PetscInt ratio;
 31:   Mat      Ii;                  /* interpolation from coarse to fine */
 32: } AppCtx;

 34: #define COARSE_LEVEL 0
 35: #define FINE_LEVEL   1

 37: /*
 38:       Mm_ratio - ration of grid lines between fine and coarse grids.
 39: */
 42: int main(int argc,char **argv)
 43: {
 45:   AppCtx         user;
 46:   PetscInt       Npx=PETSC_DECIDE,Npy=PETSC_DECIDE,Npz=PETSC_DECIDE;
 47:   PetscMPIInt    size,rank;
 48:   PetscInt       m,n,M,N,i,nrows;
 49:   PetscScalar    one = 1.0;
 50:   PetscReal      fill=2.0;
 51:   Mat            A,A_tmp,P,C,C1,C2;
 52:   PetscScalar    *array,none = -1.0,alpha;
 53:   Vec            x,v1,v2,v3,v4;
 54:   PetscReal      norm,norm_tmp,norm_tmp1,tol=1.e-12;
 55:   PetscRandom    rdm;
 56:   PetscBool      Test_MatMatMult=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_3D=PETSC_FALSE,flg;

 58:   PetscInitialize(&argc,&argv,NULL,help);
 59:   PetscOptionsGetReal(NULL,"-tol",&tol,NULL);

 61:   user.ratio     = 2;
 62:   user.coarse.mx = 2; user.coarse.my = 2; user.coarse.mz = 0;

 64:   PetscOptionsGetInt(NULL,"-Mx",&user.coarse.mx,NULL);
 65:   PetscOptionsGetInt(NULL,"-My",&user.coarse.my,NULL);
 66:   PetscOptionsGetInt(NULL,"-Mz",&user.coarse.mz,NULL);
 67:   PetscOptionsGetInt(NULL,"-ratio",&user.ratio,NULL);

 69:   if (user.coarse.mz) Test_3D = PETSC_TRUE;

 71:   user.fine.mx = user.ratio*(user.coarse.mx-1)+1;
 72:   user.fine.my = user.ratio*(user.coarse.my-1)+1;
 73:   user.fine.mz = user.ratio*(user.coarse.mz-1)+1;

 75:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 76:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 77:   PetscOptionsGetInt(NULL,"-Npx",&Npx,NULL);
 78:   PetscOptionsGetInt(NULL,"-Npy",&Npy,NULL);
 79:   PetscOptionsGetInt(NULL,"-Npz",&Npz,NULL);

 81:   /* Set up distributed array for fine grid */
 82:   if (!Test_3D) {
 83:     DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.fine.mx,
 84:                         user.fine.my,Npx,Npy,1,1,NULL,NULL,&user.fine.da);
 85:   } else {
 86:     DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,
 87:                         user.fine.mx,user.fine.my,user.fine.mz,Npx,Npy,Npz,
 88:                         1,1,NULL,NULL,NULL,&user.fine.da);
 89:   }

 91:   /* Test DMCreateMatrix()                                         */
 92:   /*------------------------------------------------------------*/
 93:   DMCreateMatrix(user.fine.da,MATAIJ,&A);
 94:   DMCreateMatrix(user.fine.da,MATBAIJ,&C);

 96:   MatConvert(C,MATAIJ,MAT_INITIAL_MATRIX,&A_tmp); /* not work for mpisbaij matrix! */
 97:   MatEqual(A,A_tmp,&flg);
 98:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"A != C");
 99:   MatDestroy(&C);
100:   MatDestroy(&A_tmp);

102:   /*------------------------------------------------------------*/

104:   MatGetLocalSize(A,&m,&n);
105:   MatGetSize(A,&M,&N);
106:   /* set val=one to A */
107:   if (size == 1) {
108:     const PetscInt *ia,*ja;
109:     MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
110:     if (flg) {
111:       MatSeqAIJGetArray(A,&array);
112:       for (i=0; i<ia[nrows]; i++) array[i] = one;
113:       MatSeqAIJRestoreArray(A,&array);
114:     }
115:     MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
116:   } else {
117:     Mat_MPIAIJ *aij = (Mat_MPIAIJ*)A->data;
118:     Mat_SeqAIJ *a   = (Mat_SeqAIJ*)(aij->A)->data, *b=(Mat_SeqAIJ*)(aij->B)->data;
119:     /* A_part */
120:     for (i=0; i<a->i[m]; i++) a->a[i] = one;
121:     /* B_part */
122:     for (i=0; i<b->i[m]; i++) b->a[i] = one;

124:   }
125:   /* MatView(A, PETSC_VIEWER_STDOUT_WORLD); */

127:   /* Set up distributed array for coarse grid */
128:   if (!Test_3D) {
129:     DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,
130:                         user.coarse.my,Npx,Npy,1,1,NULL,NULL,&user.coarse.da);
131:   } else {
132:     DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,
133:                         user.coarse.mx,user.coarse.my,user.coarse.mz,Npx,Npy,Npz,
134:                         1,1,NULL,NULL,NULL,&user.coarse.da);
135:   }

137:   /* Create interpolation between the levels */
138:   DMCreateInterpolation(user.coarse.da,user.fine.da,&P,NULL);

140:   MatGetLocalSize(P,&m,&n);
141:   MatGetSize(P,&M,&N);

143:   /* Create vectors v1 and v2 that are compatible with A */
144:   VecCreate(PETSC_COMM_WORLD,&v1);
145:   MatGetLocalSize(A,&m,NULL);
146:   VecSetSizes(v1,m,PETSC_DECIDE);
147:   VecSetFromOptions(v1);
148:   VecDuplicate(v1,&v2);
149:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
150:   PetscRandomSetFromOptions(rdm);

152:   /* Test MatMatMult(): C = A*P */
153:   /*----------------------------*/
154:   if (Test_MatMatMult) {
155:     MatDuplicate(A,MAT_COPY_VALUES,&A_tmp);
156:     MatMatMult(A_tmp,P,MAT_INITIAL_MATRIX,fill,&C);

158:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
159:     alpha=1.0;
160:     for (i=0; i<2; i++) {
161:       alpha -=0.1;
162:       MatScale(A_tmp,alpha);
163:       MatMatMult(A_tmp,P,MAT_REUSE_MATRIX,fill,&C);
164:     }

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",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:     /* Test MatDuplicate()        */
215:     /*----------------------------*/
216:     MatDuplicate(C,MAT_COPY_VALUES,&C1);
217:     MatDuplicate(C1,MAT_COPY_VALUES,&C2);
218:     MatDestroy(&C1);
219:     MatDestroy(&C2);

221:     /* Create vector x that is compatible with P */
222:     VecCreate(PETSC_COMM_WORLD,&x);
223:     MatGetLocalSize(P,&m,&n);
224:     VecSetSizes(x,n,PETSC_DECIDE);
225:     VecSetFromOptions(x);

227:     VecCreate(PETSC_COMM_WORLD,&v3);
228:     VecSetSizes(v3,n,PETSC_DECIDE);
229:     VecSetFromOptions(v3);
230:     VecDuplicate(v3,&v4);

232:     norm = 0.0;
233:     for (i=0; i<10; i++) {
234:       VecSetRandom(x,rdm);
235:       MatMult(P,x,v1);
236:       MatMult(A,v1,v2);  /* v2 = A*P*x */

238:       MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
239:       MatMult(C,x,v4);           /* v3 = C*x   */
240:       VecAXPY(v4,none,v3);
241:       VecNorm(v4,NORM_1,&norm_tmp);
242:       VecNorm(v3,NORM_1,&norm_tmp1);

244:       norm_tmp /= norm_tmp1;
245:       if (norm_tmp > norm) norm = norm_tmp;
246:     }
247:     if (norm >= tol && !rank) {
248:       PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v3 - v4|/|v3|: %G\n",norm);
249:     }
250:     MatDestroy(&C);
251:     VecDestroy(&v3);
252:     VecDestroy(&v4);
253:     VecDestroy(&x);
254:   }

256:   /* Clean up */
257:   MatDestroy(&A);
258:   PetscRandomDestroy(&rdm);
259:   VecDestroy(&v1);
260:   VecDestroy(&v2);
261:   DMDestroy(&user.fine.da);
262:   DMDestroy(&user.coarse.da);
263:   MatDestroy(&P);
264:   PetscFinalize();
265:   return 0;
266: }