Actual source code: ex5.c

  1: /*$Id: ex5.c,v 1.76 2001/08/07 21:30:37 bsmith Exp $*/

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

 13: /*
 14:   This is not a good example to understand the use of multigrid with PETSc.
 15: */
 16:  #include petscmg.h

 18: int  residual(Mat,Vec,Vec,Vec);
 19: int  gauss_seidel(void *,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,int);
 20: int  jacobi(void *,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,int);
 21: int  interpolate(Mat,Vec,Vec,Vec);
 22: int  restrct(Mat,Vec,Vec);
 23: int  Create1dLaplacian(int,Mat*);
 24: int  CalculateRhs(Vec);
 25: int  CalculateError(Vec,Vec,Vec,PetscReal*);
 26: int  CalculateSolution(int,Vec*);
 27: int  amult(Mat,Vec,Vec);

 31: int main(int Argc,char **Args)
 32: {
 33:   int         x_mesh = 15,levels = 3,cycles = 1,use_jacobi = 0;
 34:   int         i,smooths = 1,*N;
 35:   int         ierr,its;
 36:   MGType      am = MGMULTIPLICATIVE;
 37:   Mat         cmat,mat[20],fmat;
 38:   SLES        csles,sles[20],slesmg;
 39:   PetscReal   e[3]; /* l_2 error,max error, residual */
 40:   char        *shellname;
 41:   Vec         x,solution,X[20],R[20],B[20];
 42:   PetscScalar zero = 0.0;
 43:   KSP         ksp,kspmg;
 44:   PC          pcmg,pc;
 45:   PetscTruth  flg;

 47:   PetscInitialize(&Argc,&Args,(char *)0,help);

 49:   PetscOptionsGetInt(PETSC_NULL,"-x",&x_mesh,PETSC_NULL);
 50:   PetscOptionsGetInt(PETSC_NULL,"-l",&levels,PETSC_NULL);
 51:   PetscOptionsGetInt(PETSC_NULL,"-c",&cycles,PETSC_NULL);
 52:   PetscOptionsGetInt(PETSC_NULL,"-smooths",&smooths,PETSC_NULL);
 53:   PetscOptionsHasName(PETSC_NULL,"-a",&flg);
 54:   if (flg) {am = MGADDITIVE;}
 55:   PetscOptionsHasName(PETSC_NULL,"-f",&flg);
 56:   if (flg) {am = MGFULL;}
 57:   PetscOptionsHasName(PETSC_NULL,"-j",&flg);
 58:   if (flg) {use_jacobi = 1;}
 59: 
 60:   PetscMalloc(levels*sizeof(int),&N);
 61:   N[0] = x_mesh;
 62:   for (i=1; i<levels; i++) {
 63:     N[i] = N[i-1]/2;
 64:     if (N[i] < 1) {SETERRQ(1,"Too many levels");}
 65:   }

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

 69:   SLESCreate(PETSC_COMM_WORLD,&slesmg);
 70:   SLESGetPC(slesmg,&pcmg);
 71:   SLESGetKSP(slesmg,&kspmg);
 72:   SLESSetFromOptions(slesmg);
 73:   PCSetType(pcmg,PCMG);
 74:   MGSetLevels(pcmg,levels,PETSC_NULL);
 75:   MGSetType(pcmg,am);

 77:   MGGetCoarseSolve(pcmg,&csles);
 78:   SLESSetOperators(csles,cmat,cmat,DIFFERENT_NONZERO_PATTERN);
 79:   SLESGetPC(csles,&pc);
 80:   PCSetType(pc,PCLU);
 81:   SLESGetKSP(csles,&ksp);
 82:   KSPSetType(ksp,KSPPREONLY);

 84:   /* zero is finest level */
 85:   for (i=0; i<levels-1; i++) {
 86:     MGSetResidual(pcmg,levels - 1 - i,residual,(Mat)0);
 87:     MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],(void *)0,&mat[i]);
 88:     MatShellSetOperation(mat[i],MATOP_MULT,(void(*)(void))restrct);
 89:     MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void(*)(void))interpolate);
 90:     MGSetInterpolate(pcmg,levels - 1 - i,mat[i]);
 91:     MGSetRestriction(pcmg,levels - 1 - i,mat[i]);
 92:     MGSetCyclesOnLevel(pcmg,levels - 1 - i,cycles);

 94:     /* set smoother */
 95:     MGGetSmoother(pcmg,levels - 1 - i,&sles[i]);
 96:     SLESGetPC(sles[i],&pc);
 97:     PCSetType(pc,PCSHELL);
 98:     PCShellSetName(pc,"user_precond");
 99:     PCShellGetName(pc,&shellname);
100:     PetscPrintf(PETSC_COMM_WORLD,"level=%d, PCShell name is %s\n",i,shellname);

102:     /* this is a dummy! since SLES requires a matrix passed in  */
103:     SLESSetOperators(sles[i],mat[i],mat[i],DIFFERENT_NONZERO_PATTERN);
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:     PCShellSetApplyRichardson(pc,gauss_seidel,(void *)0);
110:     if (use_jacobi) {
111:       PCShellSetApplyRichardson(pc,jacobi,(void *)0);
112:     }
113:     SLESGetKSP(sles[i],&ksp);
114:     KSPSetType(ksp,KSPRICHARDSON);
115:     KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
116:     KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);

118:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
119:     X[levels - 1 - i] = x;
120:     MGSetX(pcmg,levels - 1 - i,x);
121:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
122:     B[levels -1 - i] = x;
123:     MGSetRhs(pcmg,levels - 1 - i,x);
124:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
125:     R[levels - 1 - i] = x;
126:     MGSetR(pcmg,levels - 1 - i,x);
127:   }
128:   /* create coarse level vectors */
129:   VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
130:   MGSetX(pcmg,0,x); X[0] = x;
131:   VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
132:   MGSetRhs(pcmg,0,x); B[0] = x;
133:   VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
134:   MGSetR(pcmg,0,x); R[0] = x;

136:   /* create matrix multiply for finest level */
137:   MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],(void *)0,&fmat);
138:   MatShellSetOperation(fmat,MATOP_MULT,(void(*)(void))amult);
139:   SLESSetOperators(slesmg,fmat,fmat,DIFFERENT_NONZERO_PATTERN);

141:   CalculateSolution(N[0],&solution);
142:   CalculateRhs(B[levels-1]);
143:   VecSet(&zero,X[levels-1]);

145:   if (MGCheck(pcmg)) {SETERRQ(PETSC_ERR_PLIB, "MGCheck failed");}
146: 
147:   residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
148:   CalculateError(solution,X[levels-1],R[levels-1],e);
149:   PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",e[0],e[1],e[2]);

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

156:   PetscFree(N);
157:   VecDestroy(solution);

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

176: /* --------------------------------------------------------------------- */
179: int residual(Mat mat,Vec bb,Vec xx,Vec rr)
180: {
181:   int         i,n1,ierr;
182:   PetscScalar *b,*x,*r;

185:   VecGetSize(bb,&n1);
186:   VecGetArray(bb,&b);
187:   VecGetArray(xx,&x);
188:   VecGetArray(rr,&r);
189:   n1--;
190:   r[0] = b[0] + x[1] - 2.0*x[0];
191:   r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
192:   for (i=1; i<n1; i++) {
193:     r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
194:   }
195:   VecRestoreArray(bb,&b);
196:   VecRestoreArray(xx,&x);
197:   VecRestoreArray(rr,&r);
198:   return(0);
199: }
202: int amult(Mat mat,Vec xx,Vec yy)
203: {
204:   int         i,n1,ierr;
205:   PetscScalar *y,*x;

208:   VecGetSize(xx,&n1);
209:   VecGetArray(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++) {
215:     y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
216:   }
217:   VecRestoreArray(xx,&x);
218:   VecRestoreArray(yy,&y);
219:   return(0);
220: }
221: /* --------------------------------------------------------------------- */
224: int gauss_seidel(void *ptr,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal atol,PetscReal dtol,int m)
225: {
226:   int         i,n1,ierr;
227:   PetscScalar *x,*b;

230:   VecGetSize(bb,&n1); n1--;
231:   VecGetArray(bb,&b);
232:   VecGetArray(xx,&x);
233:   while (m--) {
234:     x[0] =  .5*(x[1] + b[0]);
235:     for (i=1; i<n1; i++) {
236:       x[i] = .5*(x[i+1] + x[i-1] + b[i]);
237:     }
238:     x[n1] = .5*(x[n1-1] + b[n1]);
239:     for (i=n1-1; i>0; i--) {
240:       x[i] = .5*(x[i+1] + x[i-1] + b[i]);
241:     }
242:     x[0] =  .5*(x[1] + b[0]);
243:   }
244:   VecRestoreArray(bb,&b);
245:   VecRestoreArray(xx,&x);
246:   return(0);
247: }
248: /* --------------------------------------------------------------------- */
251: int jacobi(void *ptr,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal atol,PetscReal dtol,int m)
252: {
253:   int         i,n,n1,ierr;
254:   PetscScalar *r,*b,*x;

257:   VecGetSize(bb,&n); n1 = n - 1;
258:   VecGetArray(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++) {
265:        r[i] = .5*(x[i+1] + x[i-1] + b[i]);
266:     }
267:     r[n1] = .5*(x[n1-1] + b[n1]);
268:     for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
269:   }
270:   VecRestoreArray(bb,&b);
271:   VecRestoreArray(xx,&x);
272:   VecRestoreArray(w,&r);
273:   return(0);
274: }
275: /*
276:    We know for this application that yy  and zz are the same
277: */
278: /* --------------------------------------------------------------------- */
281: int interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
282: {
283:   int         i,n,N,i2,ierr;
284:   PetscScalar *x,*y;

287:   VecGetSize(yy,&N);
288:   VecGetArray(xx,&x);
289:   VecGetArray(yy,&y);
290:   n = N/2;
291:   for (i=0; i<n; i++) {
292:     i2 = 2*i;
293:     y[i2] +=  .5*x[i];
294:     y[i2+1] +=  x[i];
295:     y[i2+2] +=  .5*x[i];
296:   }
297:   VecRestoreArray(xx,&x);
298:   VecRestoreArray(yy,&y);
299:   return(0);
300: }
301: /* --------------------------------------------------------------------- */
304: int restrct(Mat mat,Vec rr,Vec bb)
305: {
306:   int         i,n,N,i2,ierr;
307:   PetscScalar *r,*b;

310:   VecGetSize(rr,&N);
311:   VecGetArray(rr,&r);
312:   VecGetArray(bb,&b);
313:   n = N/2;

315:   for (i=0; i<n; i++) {
316:     i2 = 2*i;
317:     b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
318:   }
319:   VecRestoreArray(rr,&r);
320:   VecRestoreArray(bb,&b);
321:   return(0);
322: }
323: /* --------------------------------------------------------------------- */
326: int Create1dLaplacian(int n,Mat *mat)
327: {
328:   PetscScalar mone = -1.0,two = 2.0;
329:   int         ierr,i,idx;

332:   MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,mat);
333: 
334:   idx= n-1;
335:   MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
336:   for (i=0; i<n-1; i++) {
337:     MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
338:     idx = i+1;
339:     MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
340:     MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
341:   }
342:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
343:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
344:   return(0);
345: }
346: /* --------------------------------------------------------------------- */
349: int CalculateRhs(Vec u)
350: {
351:   int         i,n,ierr;
352:   PetscReal   h,x = 0.0;
353:   PetscScalar uu;

356:   VecGetSize(u,&n);
357:   h = 1.0/((PetscReal)(n+1));
358:   for (i=0; i<n; i++) {
359:     x += h; uu = 2.0*h*h;
360:     VecSetValues(u,1,&i,&uu,INSERT_VALUES);
361:   }

363:   return(0);
364: }
365: /* --------------------------------------------------------------------- */
368: int CalculateSolution(int n,Vec *solution)
369: {
370:   int         i,ierr;
371:   PetscReal   h,x = 0.0;
372:   PetscScalar uu;

375:   VecCreateSeq(PETSC_COMM_SELF,n,solution);
376:   h = 1.0/((PetscReal)(n+1));
377:   for (i=0; i<n; i++) {
378:     x += h; uu = x*(1.-x);
379:     VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
380:   }
381:   return(0);
382: }
383: /* --------------------------------------------------------------------- */
386: int CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
387: {
388:   PetscScalar mone = -1.0;
389:   int    ierr;

392:   VecNorm(r,NORM_2,e+2);
393:   VecWAXPY(&mone,u,solution,r);
394:   VecNorm(r,NORM_2,e);
395:   VecNorm(r,NORM_1,e+1);
396:   return(0);
397: }