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=%" PetscInt_FMT ", 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) PCShellSetApplyRichardson(pc, jacobi_smoother);
112:     KSPSetType(ksp[i], KSPRICHARDSON);
113:     KSPSetInitialGuessNonzero(ksp[i], PETSC_TRUE);
114:     KSPSetTolerances(ksp[i], PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, smooths);

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

118:     X[levels - 1 - i] = x;
119:     if (i > 0) PCMGSetX(pcmg, levels - 1 - i, x);
120:     VecCreateSeq(PETSC_COMM_SELF, N[i], &x);

122:     B[levels - 1 - i] = x;
123:     if (i > 0) PCMGSetRhs(pcmg, levels - 1 - i, x);
124:     VecCreateSeq(PETSC_COMM_SELF, N[i], &x);

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

128:     PCMGSetR(pcmg, levels - 1 - i, x);
129:   }
130:   /* create coarse level vectors */
131:   VecCreateSeq(PETSC_COMM_SELF, N[levels - 1], &x);
132:   PCMGSetX(pcmg, 0, x);
133:   X[0] = x;
134:   VecCreateSeq(PETSC_COMM_SELF, N[levels - 1], &x);
135:   PCMGSetRhs(pcmg, 0, x);
136:   B[0] = x;

138:   /* create matrix multiply for finest level */
139:   MatCreateShell(PETSC_COMM_WORLD, N[0], N[0], N[0], N[0], NULL, &fmat);
140:   MatShellSetOperation(fmat, MATOP_MULT, (void (*)(void))amult);
141:   KSPSetOperators(kspmg, fmat, fmat);

143:   CalculateSolution(N[0], &solution);
144:   CalculateRhs(B[levels - 1]);
145:   VecSet(X[levels - 1], 0.0);

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", (double)e[0], (double)e[1], (double)e[2]);

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

157:   PetscFree(N);
158:   VecDestroy(&solution);

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

175: PetscErrorCode residual(Mat mat, Vec bb, Vec xx, Vec rr)
176: {
177:   PetscInt           i, n1;
178:   PetscScalar       *x, *r;
179:   const PetscScalar *b;

181:   VecGetSize(bb, &n1);
182:   VecGetArrayRead(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++) r[i] = b[i] + x[i + 1] + x[i - 1] - 2.0 * x[i];
189:   VecRestoreArrayRead(bb, &b);
190:   VecRestoreArray(xx, &x);
191:   VecRestoreArray(rr, &r);
192:   return 0;
193: }

195: PetscErrorCode amult(Mat mat, Vec xx, Vec yy)
196: {
197:   PetscInt           i, n1;
198:   PetscScalar       *y;
199:   const PetscScalar *x;

201:   VecGetSize(xx, &n1);
202:   VecGetArrayRead(xx, &x);
203:   VecGetArray(yy, &y);
204:   n1--;
205:   y[0]  = -x[1] + 2.0 * x[0];
206:   y[n1] = -x[n1 - 1] + 2.0 * x[n1];
207:   for (i = 1; i < n1; i++) y[i] = -x[i + 1] - x[i - 1] + 2.0 * x[i];
208:   VecRestoreArrayRead(xx, &x);
209:   VecRestoreArray(yy, &y);
210:   return 0;
211: }

213: PetscErrorCode apply_pc(PC pc, Vec bb, Vec xx)
214: {
215:   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented");
216: }

218: 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)
219: {
220:   PetscInt           i, n1;
221:   PetscScalar       *x;
222:   const PetscScalar *b;

224:   *its    = m;
225:   *reason = PCRICHARDSON_CONVERGED_ITS;
226:   VecGetSize(bb, &n1);
227:   n1--;
228:   VecGetArrayRead(bb, &b);
229:   VecGetArray(xx, &x);
230:   while (m--) {
231:     x[0] = .5 * (x[1] + b[0]);
232:     for (i = 1; i < n1; i++) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
233:     x[n1] = .5 * (x[n1 - 1] + b[n1]);
234:     for (i = n1 - 1; i > 0; i--) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
235:     x[0] = .5 * (x[1] + b[0]);
236:   }
237:   VecRestoreArrayRead(bb, &b);
238:   VecRestoreArray(xx, &x);
239:   return 0;
240: }

242: 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)
243: {
244:   PetscInt           i, n, n1;
245:   PetscScalar       *r, *x;
246:   const PetscScalar *b;

248:   *its    = m;
249:   *reason = PCRICHARDSON_CONVERGED_ITS;
250:   VecGetSize(bb, &n);
251:   n1 = n - 1;
252:   VecGetArrayRead(bb, &b);
253:   VecGetArray(xx, &x);
254:   VecGetArray(w, &r);

256:   while (m--) {
257:     r[0] = .5 * (x[1] + b[0]);
258:     for (i = 1; i < n1; i++) r[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
259:     r[n1] = .5 * (x[n1 - 1] + b[n1]);
260:     for (i = 0; i < n; i++) x[i] = (2.0 * r[i] + x[i]) / 3.0;
261:   }
262:   VecRestoreArrayRead(bb, &b);
263:   VecRestoreArray(xx, &x);
264:   VecRestoreArray(w, &r);
265:   return 0;
266: }
267: /*
268:    We know for this application that yy  and zz are the same
269: */

271: PetscErrorCode interpolate(Mat mat, Vec xx, Vec yy, Vec zz)
272: {
273:   PetscInt           i, n, N, i2;
274:   PetscScalar       *y;
275:   const PetscScalar *x;

277:   VecGetSize(yy, &N);
278:   VecGetArrayRead(xx, &x);
279:   VecGetArray(yy, &y);
280:   n = N / 2;
281:   for (i = 0; i < n; i++) {
282:     i2 = 2 * i;
283:     y[i2] += .5 * x[i];
284:     y[i2 + 1] += x[i];
285:     y[i2 + 2] += .5 * x[i];
286:   }
287:   VecRestoreArrayRead(xx, &x);
288:   VecRestoreArray(yy, &y);
289:   return 0;
290: }

292: PetscErrorCode restrct(Mat mat, Vec rr, Vec bb)
293: {
294:   PetscInt           i, n, N, i2;
295:   PetscScalar       *b;
296:   const PetscScalar *r;

298:   VecGetSize(rr, &N);
299:   VecGetArrayRead(rr, &r);
300:   VecGetArray(bb, &b);
301:   n = N / 2;

303:   for (i = 0; i < n; i++) {
304:     i2   = 2 * i;
305:     b[i] = (r[i2] + 2.0 * r[i2 + 1] + r[i2 + 2]);
306:   }
307:   VecRestoreArrayRead(rr, &r);
308:   VecRestoreArray(bb, &b);
309:   return 0;
310: }

312: PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat)
313: {
314:   PetscScalar mone = -1.0, two = 2.0;
315:   PetscInt    i, idx;

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

319:   idx = n - 1;
320:   MatSetValues(*mat, 1, &idx, 1, &idx, &two, INSERT_VALUES);
321:   for (i = 0; i < n - 1; i++) {
322:     MatSetValues(*mat, 1, &i, 1, &i, &two, INSERT_VALUES);
323:     idx = i + 1;
324:     MatSetValues(*mat, 1, &idx, 1, &i, &mone, INSERT_VALUES);
325:     MatSetValues(*mat, 1, &i, 1, &idx, &mone, INSERT_VALUES);
326:   }
327:   MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
328:   MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
329:   return 0;
330: }

332: PetscErrorCode CalculateRhs(Vec u)
333: {
334:   PetscInt    i, n;
335:   PetscReal   h;
336:   PetscScalar uu;

338:   VecGetSize(u, &n);
339:   h = 1.0 / ((PetscReal)(n + 1));
340:   for (i = 0; i < n; i++) {
341:     uu = 2.0 * h * h;
342:     VecSetValues(u, 1, &i, &uu, INSERT_VALUES);
343:   }
344:   return 0;
345: }

347: PetscErrorCode CalculateSolution(PetscInt n, Vec *solution)
348: {
349:   PetscInt    i;
350:   PetscReal   h, x = 0.0;
351:   PetscScalar uu;

353:   VecCreateSeq(PETSC_COMM_SELF, n, solution);
354:   h = 1.0 / ((PetscReal)(n + 1));
355:   for (i = 0; i < n; i++) {
356:     x += h;
357:     uu = x * (1. - x);
358:     VecSetValues(*solution, 1, &i, &uu, INSERT_VALUES);
359:   }
360:   return 0;
361: }

363: PetscErrorCode CalculateError(Vec solution, Vec u, Vec r, PetscReal *e)
364: {
365:   VecNorm(r, NORM_2, e + 2);
366:   VecWAXPY(r, -1.0, u, solution);
367:   VecNorm(r, NORM_2, e);
368:   VecNorm(r, NORM_1, e + 1);
369:   return 0;
370: }

372: /*TEST

374:    test:

376: TEST*/