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: */
 15:  #include petscmg.h

 17: PetscErrorCode  residual(Mat,Vec,Vec,Vec);
 18: PetscErrorCode  gauss_seidel(PC,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,PetscInt,PetscTruth,PetscInt*,PCRichardsonConvergedReason*);
 19: PetscErrorCode  jacobi(PC,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,PetscInt,PetscTruth,PetscInt*,PCRichardsonConvergedReason*);
 20: PetscErrorCode  interpolate(Mat,Vec,Vec,Vec);
 21: PetscErrorCode  restrct(Mat,Vec,Vec);
 22: PetscErrorCode  Create1dLaplacian(PetscInt,Mat*);
 23: PetscErrorCode  CalculateRhs(Vec);
 24: PetscErrorCode  CalculateError(Vec,Vec,Vec,PetscReal*);
 25: PetscErrorCode  CalculateSolution(PetscInt,Vec*);
 26: PetscErrorCode  amult(Mat,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:   PetscErrorCode  ierr;
 35:   PCMGType        am = PC_MG_MULTIPLICATIVE;
 36:   Mat             cmat,mat[20],fmat;
 37:   KSP             cksp,ksp[20],kspmg;
 38:   PetscReal       e[3]; /* l_2 error,max error, residual */
 39:   char            *shellname;
 40:   Vec             x,solution,X[20],R[20],B[20];
 41:   PC              pcmg,pc;
 42:   PetscTruth      flg;

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

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

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

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

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

 79:   /* zero is finest level */
 80:   for (i=0; i<levels-1; i++) {
 81:     PCMGSetResidual(pcmg,levels - 1 - i,residual,(Mat)0);
 82:     MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],(void*)0,&mat[i]);
 83:     MatShellSetOperation(mat[i],MATOP_MULT,(void(*)(void))restrct);
 84:     MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void(*)(void))interpolate);
 85:     PCMGSetInterpolation(pcmg,levels - 1 - i,mat[i]);
 86:     PCMGSetRestriction(pcmg,levels - 1 - i,mat[i]);
 87:     PCMGSetCyclesOnLevel(pcmg,levels - 1 - i,cycles);

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

 97:     /* this is a dummy! since KSP requires a matrix passed in  */
 98:     KSPSetOperators(ksp[i],mat[i],mat[i],DIFFERENT_NONZERO_PATTERN);
 99:     /* 
100:         We override the matrix passed in by forcing it to use Richardson with 
101:         a user provided application. This is non-standard and this practice
102:         should be avoided.
103:     */
104:     PCShellSetApplyRichardson(pc,gauss_seidel);
105:     if (use_jacobi) {
106:       PCShellSetApplyRichardson(pc,jacobi);
107:     }
108:     KSPSetType(ksp[i],KSPRICHARDSON);
109:     KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);
110:     KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);

112:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
113:     X[levels - 1 - i] = x;
114:     if (i > 0) {
115:       PCMGSetX(pcmg,levels - 1 - i,x);
116:     }
117:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
118:     B[levels -1 - i] = x;
119:     if (i > 0) {
120:       PCMGSetRhs(pcmg,levels - 1 - i,x);
121:     }
122:     VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
123:     R[levels - 1 - i] = x;
124:     PCMGSetR(pcmg,levels - 1 - i,x);
125:   }
126:   /* create coarse level vectors */
127:   VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
128:   PCMGSetX(pcmg,0,x); X[0] = x;
129:   VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
130:   PCMGSetRhs(pcmg,0,x); B[0] = x;

132:   /* create matrix multiply for finest level */
133:   MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],(void*)0,&fmat);
134:   MatShellSetOperation(fmat,MATOP_MULT,(void(*)(void))amult);
135:   KSPSetOperators(kspmg,fmat,fmat,DIFFERENT_NONZERO_PATTERN);

137:   CalculateSolution(N[0],&solution);
138:   CalculateRhs(B[levels-1]);
139:   VecSet(X[levels-1],0.0);

141:   residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
142:   CalculateError(solution,X[levels-1],R[levels-1],e);
143:   PetscPrintf(PETSC_COMM_SELF,"l_2 error %G max error %G resi %G\n",e[0],e[1],e[2]);

145:   KSPSolve(kspmg,B[levels-1],X[levels-1]);
146:   KSPGetIterationNumber(kspmg,&its);
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,"its %D l_2 error %G max error %G resi %G\n",its,e[0],e[1],e[2]);

151:   PetscFree(N);
152:   VecDestroy(solution);

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

171: /* --------------------------------------------------------------------- */
174: PetscErrorCode residual(Mat mat,Vec bb,Vec xx,Vec rr)
175: {
176:   PetscInt       i,n1;
178:   PetscScalar    *b,*x,*r;

181:   VecGetSize(bb,&n1);
182:   VecGetArray(bb,&b);
183:   VecGetArray(xx,&x);
184:   VecGetArray(rr,&r);
185:   n1--;
186:   r[0] = b[0] + x[1] - 2.0*x[0];
187:   r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
188:   for (i=1; i<n1; i++) {
189:     r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
190:   }
191:   VecRestoreArray(bb,&b);
192:   VecRestoreArray(xx,&x);
193:   VecRestoreArray(rr,&r);
194:   return(0);
195: }
198: PetscErrorCode amult(Mat mat,Vec xx,Vec yy)
199: {
200:   PetscInt       i,n1;
202:   PetscScalar    *y,*x;

205:   VecGetSize(xx,&n1);
206:   VecGetArray(xx,&x);
207:   VecGetArray(yy,&y);
208:   n1--;
209:   y[0] =  -x[1] + 2.0*x[0];
210:   y[n1] = -x[n1-1] + 2.0*x[n1];
211:   for (i=1; i<n1; i++) {
212:     y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
213:   }
214:   VecRestoreArray(xx,&x);
215:   VecRestoreArray(yy,&y);
216:   return(0);
217: }
218: /* --------------------------------------------------------------------- */
221: PetscErrorCode gauss_seidel(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscTruth guesszero,PetscInt *its,PCRichardsonConvergedReason *reason)
222: {
223:   PetscInt       i,n1;
225:   PetscScalar    *x,*b;

228:   *its    = m;
229:   *reason = PCRICHARDSON_CONVERGED_ITS;
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: PetscErrorCode jacobi(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscTruth guesszero,PetscInt *its,PCRichardsonConvergedReason *reason)
252: {
253:   PetscInt       i,n,n1;
255:   PetscScalar    *r,*b,*x;

258:   *its = m;
259:   *reason = PCRICHARDSON_CONVERGED_ITS;
260:   VecGetSize(bb,&n); n1 = n - 1;
261:   VecGetArray(bb,&b);
262:   VecGetArray(xx,&x);
263:   VecGetArray(w,&r);

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

291:   VecGetSize(yy,&N);
292:   VecGetArray(xx,&x);
293:   VecGetArray(yy,&y);
294:   n = N/2;
295:   for (i=0; i<n; i++) {
296:     i2 = 2*i;
297:     y[i2] +=  .5*x[i];
298:     y[i2+1] +=  x[i];
299:     y[i2+2] +=  .5*x[i];
300:   }
301:   VecRestoreArray(xx,&x);
302:   VecRestoreArray(yy,&y);
303:   return(0);
304: }
305: /* --------------------------------------------------------------------- */
308: PetscErrorCode restrct(Mat mat,Vec rr,Vec bb)
309: {
310:   PetscInt       i,n,N,i2;
312:   PetscScalar    *r,*b;

315:   VecGetSize(rr,&N);
316:   VecGetArray(rr,&r);
317:   VecGetArray(bb,&b);
318:   n = N/2;

320:   for (i=0; i<n; i++) {
321:     i2 = 2*i;
322:     b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
323:   }
324:   VecRestoreArray(rr,&r);
325:   VecRestoreArray(bb,&b);
326:   return(0);
327: }
328: /* --------------------------------------------------------------------- */
331: PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
332: {
333:   PetscScalar    mone = -1.0,two = 2.0;
334:   PetscInt       i,idx;

338:   MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,mat);
339: 
340:   idx= n-1;
341:   MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
342:   for (i=0; i<n-1; i++) {
343:     MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
344:     idx = i+1;
345:     MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
346:     MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
347:   }
348:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
349:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
350:   return(0);
351: }
352: /* --------------------------------------------------------------------- */
355: PetscErrorCode CalculateRhs(Vec u)
356: {
358:   PetscInt       i,n;
359:   PetscReal      h,x = 0.0;
360:   PetscScalar    uu;

363:   VecGetSize(u,&n);
364:   h = 1.0/((PetscReal)(n+1));
365:   for (i=0; i<n; i++) {
366:     x += h; uu = 2.0*h*h;
367:     VecSetValues(u,1,&i,&uu,INSERT_VALUES);
368:   }

370:   return(0);
371: }
372: /* --------------------------------------------------------------------- */
375: PetscErrorCode CalculateSolution(PetscInt n,Vec *solution)
376: {
378:   PetscInt       i;
379:   PetscReal      h,x = 0.0;
380:   PetscScalar    uu;

383:   VecCreateSeq(PETSC_COMM_SELF,n,solution);
384:   h = 1.0/((PetscReal)(n+1));
385:   for (i=0; i<n; i++) {
386:     x += h; uu = x*(1.-x);
387:     VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
388:   }
389:   return(0);
390: }
391: /* --------------------------------------------------------------------- */
394: PetscErrorCode CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
395: {

399:   VecNorm(r,NORM_2,e+2);
400:   VecWAXPY(r,-1.0,u,solution);
401:   VecNorm(r,NORM_2,e);
402:   VecNorm(r,NORM_1,e+1);
403:   return(0);
404: }