Actual source code: ex5.c


  2: static char help[] = "Tests the multigrid code.  The input parameters are:\n\
  3:   -x N              Use a mesh in the x direction of N.  \n\
  4:   -c N              Use N V-cycles.  \n\
  5:   -l N              Use N Levels.  \n\
  6:   -smooths N        Use N pre smooths and N post smooths.  \n\
  7:   -j                Use Jacobi smoother.  \n\
  8:   -a use additive multigrid \n\
  9:   -f use full multigrid (preconditioner variant) \n\
 10: This example also demonstrates matrix-free methods\n\n";

 12: /*
 13:   This is not a good example to understand the use of multigrid with PETSc.
 14: */

 16: #include <petscksp.h>

 18: PetscErrorCode  residual(Mat,Vec,Vec,Vec);
 19: PetscErrorCode  gauss_seidel(PC,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,PetscInt,PetscBool,PetscInt*,PCRichardsonConvergedReason*);
 20: PetscErrorCode  jacobi_smoother(PC,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,PetscInt,PetscBool,PetscInt*,PCRichardsonConvergedReason*);
 21: PetscErrorCode  interpolate(Mat,Vec,Vec,Vec);
 22: PetscErrorCode  restrct(Mat,Vec,Vec);
 23: PetscErrorCode  Create1dLaplacian(PetscInt,Mat*);
 24: PetscErrorCode  CalculateRhs(Vec);
 25: PetscErrorCode  CalculateError(Vec,Vec,Vec,PetscReal*);
 26: PetscErrorCode  CalculateSolution(PetscInt,Vec*);
 27: PetscErrorCode  amult(Mat,Vec,Vec);
 28: PetscErrorCode  apply_pc(PC,Vec,Vec);

 30: int main(int Argc,char **Args)
 31: {
 32:   PetscInt       x_mesh = 15,levels = 3,cycles = 1,use_jacobi = 0;
 33:   PetscInt       i,smooths = 1,*N,its;
 34:   PCMGType       am = PC_MG_MULTIPLICATIVE;
 35:   Mat            cmat,mat[20],fmat;
 36:   KSP            cksp,ksp[20],kspmg;
 37:   PetscReal      e[3];  /* l_2 error,max error, residual */
 38:   const char     *shellname;
 39:   Vec            x,solution,X[20],R[20],B[20];
 40:   PC             pcmg,pc;
 41:   PetscBool      flg;

 43:   PetscInitialize(&Argc,&Args,(char*)0,help);
 44:   PetscOptionsGetInt(NULL,NULL,"-x",&x_mesh,NULL);
 45:   PetscOptionsGetInt(NULL,NULL,"-l",&levels,NULL);
 46:   PetscOptionsGetInt(NULL,NULL,"-c",&cycles,NULL);
 47:   PetscOptionsGetInt(NULL,NULL,"-smooths",&smooths,NULL);
 48:   PetscOptionsHasName(NULL,NULL,"-a",&flg);

 50:   if (flg) am = PC_MG_ADDITIVE;
 51:   PetscOptionsHasName(NULL,NULL,"-f",&flg);
 52:   if (flg) am = PC_MG_FULL;
 53:   PetscOptionsHasName(NULL,NULL,"-j",&flg);
 54:   if (flg) use_jacobi = 1;

 56:   PetscMalloc1(levels,&N);
 57:   N[0] = x_mesh;
 58:   for (i=1; i<levels; i++) {
 59:     N[i] = N[i-1]/2;
 61:   }

 63:   Create1dLaplacian(N[levels-1],&cmat);

 65:   KSPCreate(PETSC_COMM_WORLD,&kspmg);
 66:   KSPGetPC(kspmg,&pcmg);
 67:   KSPSetFromOptions(kspmg);
 68:   PCSetType(pcmg,PCMG);
 69:   PCMGSetLevels(pcmg,levels,NULL);
 70:   PCMGSetType(pcmg,am);

 72:   PCMGGetCoarseSolve(pcmg,&cksp);
 73:   KSPSetOperators(cksp,cmat,cmat);
 74:   KSPGetPC(cksp,&pc);
 75:   PCSetType(pc,PCLU);
 76:   KSPSetType(cksp,KSPPREONLY);

 78:   /* zero is finest level */
 79:   for (i=0; i<levels-1; i++) {
 80:     Mat dummy;

 82:     PCMGSetResidual(pcmg,levels - 1 - i,residual,NULL);
 83:     MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],NULL,&mat[i]);
 84:     MatShellSetOperation(mat[i],MATOP_MULT,(void (*)(void))restrct);
 85:     MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void (*)(void))interpolate);
 86:     PCMGSetInterpolation(pcmg,levels - 1 - i,mat[i]);
 87:     PCMGSetRestriction(pcmg,levels - 1 - i,mat[i]);
 88:     PCMGSetCycleTypeOnLevel(pcmg,levels - 1 - i,(PCMGCycleType)cycles);

 90:     /* set smoother */
 91:     PCMGGetSmoother(pcmg,levels - 1 - i,&ksp[i]);
 92:     KSPGetPC(ksp[i],&pc);
 93:     PCSetType(pc,PCSHELL);
 94:     PCShellSetName(pc,"user_precond");
 95:     PCShellGetName(pc,&shellname);
 96:     PetscPrintf(PETSC_COMM_WORLD,"level=%D, PCShell name is %s\n",i,shellname);

 98:     /* this is not used unless different options are passed to the solver */
 99:     MatCreateShell(PETSC_COMM_WORLD,N[i],N[i],N[i],N[i],NULL,&dummy);
100:     MatShellSetOperation(dummy,MATOP_MULT,(void (*)(void))amult);
101:     KSPSetOperators(ksp[i],dummy,dummy);
102:     MatDestroy(&dummy);

104:     /*
105:         We override the matrix passed in by forcing it to use Richardson with
106:         a user provided application. This is non-standard and this practice
107:         should be avoided.
108:     */
109:     PCShellSetApply(pc,apply_pc);
110:     PCShellSetApplyRichardson(pc,gauss_seidel);
111:     if (use_jacobi) {
112:       PCShellSetApplyRichardson(pc,jacobi_smoother);
113:     }
114:     KSPSetType(ksp[i],KSPRICHARDSON);
115:     KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);
116:     KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);

118:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);

120:     X[levels - 1 - i] = x;
121:     if (i > 0) {
122:       PCMGSetX(pcmg,levels - 1 - i,x);
123:     }
124:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);

126:     B[levels -1 - i] = x;
127:     if (i > 0) {
128:       PCMGSetRhs(pcmg,levels - 1 - i,x);
129:     }
130:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);

132:     R[levels - 1 - i] = x;

134:     PCMGSetR(pcmg,levels - 1 - i,x);
135:   }
136:   /* create coarse level vectors */
137:   VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
138:   PCMGSetX(pcmg,0,x); X[0] = x;
139:   VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
140:   PCMGSetRhs(pcmg,0,x); B[0] = x;

142:   /* create matrix multiply for finest level */
143:   MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],NULL,&fmat);
144:   MatShellSetOperation(fmat,MATOP_MULT,(void (*)(void))amult);
145:   KSPSetOperators(kspmg,fmat,fmat);

147:   CalculateSolution(N[0],&solution);
148:   CalculateRhs(B[levels-1]);
149:   VecSet(X[levels-1],0.0);

151:   residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
152:   CalculateError(solution,X[levels-1],R[levels-1],e);
153:   PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",(double)e[0],(double)e[1],(double)e[2]);

155:   KSPSolve(kspmg,B[levels-1],X[levels-1]);
156:   KSPGetIterationNumber(kspmg,&its);
157:   residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
158:   CalculateError(solution,X[levels-1],R[levels-1],e);
159:   PetscPrintf(PETSC_COMM_SELF,"its %D l_2 error %g max error %g resi %g\n",its,(double)e[0],(double)e[1],(double)e[2]);

161:   PetscFree(N);
162:   VecDestroy(&solution);

164:   /* note we have to keep a list of all vectors allocated, this is
165:      not ideal, but putting it in MGDestroy is not so good either*/
166:   for (i=0; i<levels; i++) {
167:     VecDestroy(&X[i]);
168:     VecDestroy(&B[i]);
169:     if (i) VecDestroy(&R[i]);
170:   }
171:   for (i=0; i<levels-1; i++) {
172:     MatDestroy(&mat[i]);
173:   }
174:   MatDestroy(&cmat);
175:   MatDestroy(&fmat);
176:   KSPDestroy(&kspmg);
177:   PetscFinalize();
178:   return 0;
179: }

181: /* --------------------------------------------------------------------- */
182: PetscErrorCode residual(Mat mat,Vec bb,Vec xx,Vec rr)
183: {
184:   PetscInt          i,n1;
185:   PetscScalar       *x,*r;
186:   const PetscScalar *b;

188:   VecGetSize(bb,&n1);
189:   VecGetArrayRead(bb,&b);
190:   VecGetArray(xx,&x);
191:   VecGetArray(rr,&r);
192:   n1--;
193:   r[0]  = b[0] + x[1] - 2.0*x[0];
194:   r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
195:   for (i=1; i<n1; i++) r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
196:   VecRestoreArrayRead(bb,&b);
197:   VecRestoreArray(xx,&x);
198:   VecRestoreArray(rr,&r);
199:   return 0;
200: }

202: PetscErrorCode amult(Mat mat,Vec xx,Vec yy)
203: {
204:   PetscInt          i,n1;
205:   PetscScalar       *y;
206:   const PetscScalar *x;

208:   VecGetSize(xx,&n1);
209:   VecGetArrayRead(xx,&x);
210:   VecGetArray(yy,&y);
211:   n1--;
212:   y[0] =  -x[1] + 2.0*x[0];
213:   y[n1] = -x[n1-1] + 2.0*x[n1];
214:   for (i=1; i<n1; i++) y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
215:   VecRestoreArrayRead(xx,&x);
216:   VecRestoreArray(yy,&y);
217:   return 0;
218: }

220: /* --------------------------------------------------------------------- */
221: PetscErrorCode apply_pc(PC pc,Vec bb,Vec xx)
222: {
223:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not implemented");
224: }

226: PetscErrorCode gauss_seidel(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscBool guesszero,PetscInt *its,PCRichardsonConvergedReason *reason)
227: {
228:   PetscInt          i,n1;
229:   PetscScalar       *x;
230:   const PetscScalar *b;

232:   *its    = m;
233:   *reason = PCRICHARDSON_CONVERGED_ITS;
234:   VecGetSize(bb,&n1); n1--;
235:   VecGetArrayRead(bb,&b);
236:   VecGetArray(xx,&x);
237:   while (m--) {
238:     x[0] =  .5*(x[1] + b[0]);
239:     for (i=1; i<n1; i++) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
240:     x[n1] = .5*(x[n1-1] + b[n1]);
241:     for (i=n1-1; i>0; i--) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
242:     x[0] =  .5*(x[1] + b[0]);
243:   }
244:   VecRestoreArrayRead(bb,&b);
245:   VecRestoreArray(xx,&x);
246:   return 0;
247: }
248: /* --------------------------------------------------------------------- */
249: PetscErrorCode jacobi_smoother(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscBool guesszero,PetscInt *its,PCRichardsonConvergedReason *reason)
250: {
251:   PetscInt          i,n,n1;
252:   PetscScalar       *r,*x;
253:   const PetscScalar *b;

255:   *its    = m;
256:   *reason = PCRICHARDSON_CONVERGED_ITS;
257:   VecGetSize(bb,&n); n1 = n - 1;
258:   VecGetArrayRead(bb,&b);
259:   VecGetArray(xx,&x);
260:   VecGetArray(w,&r);

262:   while (m--) {
263:     r[0] = .5*(x[1] + b[0]);
264:     for (i=1; i<n1; i++) r[i] = .5*(x[i+1] + x[i-1] + b[i]);
265:     r[n1] = .5*(x[n1-1] + b[n1]);
266:     for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
267:   }
268:   VecRestoreArrayRead(bb,&b);
269:   VecRestoreArray(xx,&x);
270:   VecRestoreArray(w,&r);
271:   return 0;
272: }
273: /*
274:    We know for this application that yy  and zz are the same
275: */
276: /* --------------------------------------------------------------------- */
277: PetscErrorCode interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
278: {
279:   PetscInt          i,n,N,i2;
280:   PetscScalar       *y;
281:   const PetscScalar *x;

283:   VecGetSize(yy,&N);
284:   VecGetArrayRead(xx,&x);
285:   VecGetArray(yy,&y);
286:   n    = N/2;
287:   for (i=0; i<n; i++) {
288:     i2       = 2*i;
289:     y[i2]   += .5*x[i];
290:     y[i2+1] +=    x[i];
291:     y[i2+2] += .5*x[i];
292:   }
293:   VecRestoreArrayRead(xx,&x);
294:   VecRestoreArray(yy,&y);
295:   return 0;
296: }
297: /* --------------------------------------------------------------------- */
298: PetscErrorCode restrct(Mat mat,Vec rr,Vec bb)
299: {
300:   PetscInt          i,n,N,i2;
301:   PetscScalar       *b;
302:   const PetscScalar *r;

304:   VecGetSize(rr,&N);
305:   VecGetArrayRead(rr,&r);
306:   VecGetArray(bb,&b);
307:   n    = N/2;

309:   for (i=0; i<n; i++) {
310:     i2   = 2*i;
311:     b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
312:   }
313:   VecRestoreArrayRead(rr,&r);
314:   VecRestoreArray(bb,&b);
315:   return 0;
316: }
317: /* --------------------------------------------------------------------- */
318: PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
319: {
320:   PetscScalar    mone = -1.0,two = 2.0;
321:   PetscInt       i,idx;

323:   MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,mat);

325:   idx  = n-1;
326:   MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
327:   for (i=0; i<n-1; i++) {
328:     MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
329:     idx  = i+1;
330:     MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
331:     MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
332:   }
333:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
334:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
335:   return 0;
336: }
337: /* --------------------------------------------------------------------- */
338: PetscErrorCode CalculateRhs(Vec u)
339: {
340:   PetscInt    i,n;
341:   PetscReal   h;
342:   PetscScalar uu;

344:   VecGetSize(u,&n);
345:   h    = 1.0/((PetscReal)(n+1));
346:   for (i=0; i<n; i++) {
347:     uu = 2.0*h*h;
348:     VecSetValues(u,1,&i,&uu,INSERT_VALUES);
349:   }
350:   return 0;
351: }
352: /* --------------------------------------------------------------------- */
353: PetscErrorCode CalculateSolution(PetscInt n,Vec *solution)
354: {
355:   PetscInt       i;
356:   PetscReal      h,x = 0.0;
357:   PetscScalar    uu;

359:   VecCreateSeq(PETSC_COMM_SELF,n,solution);
360:   h    = 1.0/((PetscReal)(n+1));
361:   for (i=0; i<n; i++) {
362:     x   += h; uu = x*(1.-x);
363:     VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
364:   }
365:   return 0;
366: }
367: /* --------------------------------------------------------------------- */
368: PetscErrorCode CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
369: {
370:   VecNorm(r,NORM_2,e+2);
371:   VecWAXPY(r,-1.0,u,solution);
372:   VecNorm(r,NORM_2,e);
373:   VecNorm(r,NORM_1,e+1);
374:   return 0;
375: }

377: /*TEST

379:    test:

381: TEST*/