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:   PetscCall(PetscInitialize(&Argc, &Args, (char *)0, help));
 44:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-x", &x_mesh, NULL));
 45:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-l", &levels, NULL));
 46:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-c", &cycles, NULL));
 47:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-smooths", &smooths, NULL));
 48:   PetscCall(PetscOptionsHasName(NULL, NULL, "-a", &flg));

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

 56:   PetscCall(PetscMalloc1(levels, &N));
 57:   N[0] = x_mesh;
 58:   for (i = 1; i < levels; i++) {
 59:     N[i] = N[i - 1] / 2;
 60:     PetscCheck(N[i] >= 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "Too many levels or N is not large enough");
 61:   }

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

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

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

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

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

 90:     /* set smoother */
 91:     PetscCall(PCMGGetSmoother(pcmg, levels - 1 - i, &ksp[i]));
 92:     PetscCall(KSPGetPC(ksp[i], &pc));
 93:     PetscCall(PCSetType(pc, PCSHELL));
 94:     PetscCall(PCShellSetName(pc, "user_precond"));
 95:     PetscCall(PCShellGetName(pc, &shellname));
 96:     PetscCall(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:     PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[i], N[i], N[i], N[i], NULL, &dummy));
100:     PetscCall(MatShellSetOperation(dummy, MATOP_MULT, (void (*)(void))amult));
101:     PetscCall(KSPSetOperators(ksp[i], dummy, dummy));
102:     PetscCall(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:     PetscCall(PCShellSetApply(pc, apply_pc));
110:     PetscCall(PCShellSetApplyRichardson(pc, gauss_seidel));
111:     if (use_jacobi) PetscCall(PCShellSetApplyRichardson(pc, jacobi_smoother));
112:     PetscCall(KSPSetType(ksp[i], KSPRICHARDSON));
113:     PetscCall(KSPSetInitialGuessNonzero(ksp[i], PETSC_TRUE));
114:     PetscCall(KSPSetTolerances(ksp[i], PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, smooths));

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

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

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

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

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

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

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

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

151:   PetscCall(KSPSolve(kspmg, B[levels - 1], X[levels - 1]));
152:   PetscCall(KSPGetIterationNumber(kspmg, &its));
153:   PetscCall(residual((Mat)0, B[levels - 1], X[levels - 1], R[levels - 1]));
154:   PetscCall(CalculateError(solution, X[levels - 1], R[levels - 1], e));
155:   PetscCall(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:   PetscCall(PetscFree(N));
158:   PetscCall(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:     PetscCall(VecDestroy(&X[i]));
164:     PetscCall(VecDestroy(&B[i]));
165:     if (i) PetscCall(VecDestroy(&R[i]));
166:   }
167:   for (i = 0; i < levels - 1; i++) PetscCall(MatDestroy(&mat[i]));
168:   PetscCall(MatDestroy(&cmat));
169:   PetscCall(MatDestroy(&fmat));
170:   PetscCall(KSPDestroy(&kspmg));
171:   PetscCall(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:   PetscFunctionBegin;
182:   PetscCall(VecGetSize(bb, &n1));
183:   PetscCall(VecGetArrayRead(bb, &b));
184:   PetscCall(VecGetArray(xx, &x));
185:   PetscCall(VecGetArray(rr, &r));
186:   n1--;
187:   r[0]  = b[0] + x[1] - 2.0 * x[0];
188:   r[n1] = b[n1] + x[n1 - 1] - 2.0 * x[n1];
189:   for (i = 1; i < n1; i++) r[i] = b[i] + x[i + 1] + x[i - 1] - 2.0 * x[i];
190:   PetscCall(VecRestoreArrayRead(bb, &b));
191:   PetscCall(VecRestoreArray(xx, &x));
192:   PetscCall(VecRestoreArray(rr, &r));
193:   PetscFunctionReturn(PETSC_SUCCESS);
194: }

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

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

215: PetscErrorCode apply_pc(PC pc, Vec bb, Vec xx)
216: {
217:   PetscFunctionBegin;
218:   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented");
219: }

221: 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)
222: {
223:   PetscInt           i, n1;
224:   PetscScalar       *x;
225:   const PetscScalar *b;

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

246: 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)
247: {
248:   PetscInt           i, n, n1;
249:   PetscScalar       *r, *x;
250:   const PetscScalar *b;

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

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

276: PetscErrorCode interpolate(Mat mat, Vec xx, Vec yy, Vec zz)
277: {
278:   PetscInt           i, n, N, i2;
279:   PetscScalar       *y;
280:   const PetscScalar *x;

282:   PetscFunctionBegin;
283:   PetscCall(VecGetSize(yy, &N));
284:   PetscCall(VecGetArrayRead(xx, &x));
285:   PetscCall(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:   PetscCall(VecRestoreArrayRead(xx, &x));
294:   PetscCall(VecRestoreArray(yy, &y));
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: PetscErrorCode restrct(Mat mat, Vec rr, Vec bb)
299: {
300:   PetscInt           i, n, N, i2;
301:   PetscScalar       *b;
302:   const PetscScalar *r;

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

310:   for (i = 0; i < n; i++) {
311:     i2   = 2 * i;
312:     b[i] = (r[i2] + 2.0 * r[i2 + 1] + r[i2 + 2]);
313:   }
314:   PetscCall(VecRestoreArrayRead(rr, &r));
315:   PetscCall(VecRestoreArray(bb, &b));
316:   PetscFunctionReturn(PETSC_SUCCESS);
317: }

319: PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat)
320: {
321:   PetscScalar mone = -1.0, two = 2.0;
322:   PetscInt    i, idx;

324:   PetscFunctionBegin;
325:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat));

327:   idx = n - 1;
328:   PetscCall(MatSetValues(*mat, 1, &idx, 1, &idx, &two, INSERT_VALUES));
329:   for (i = 0; i < n - 1; i++) {
330:     PetscCall(MatSetValues(*mat, 1, &i, 1, &i, &two, INSERT_VALUES));
331:     idx = i + 1;
332:     PetscCall(MatSetValues(*mat, 1, &idx, 1, &i, &mone, INSERT_VALUES));
333:     PetscCall(MatSetValues(*mat, 1, &i, 1, &idx, &mone, INSERT_VALUES));
334:   }
335:   PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
336:   PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: PetscErrorCode CalculateRhs(Vec u)
341: {
342:   PetscInt    i, n;
343:   PetscReal   h;
344:   PetscScalar uu;

346:   PetscFunctionBegin;
347:   PetscCall(VecGetSize(u, &n));
348:   h = 1.0 / ((PetscReal)(n + 1));
349:   for (i = 0; i < n; i++) {
350:     uu = 2.0 * h * h;
351:     PetscCall(VecSetValues(u, 1, &i, &uu, INSERT_VALUES));
352:   }
353:   PetscFunctionReturn(PETSC_SUCCESS);
354: }

356: PetscErrorCode CalculateSolution(PetscInt n, Vec *solution)
357: {
358:   PetscInt    i;
359:   PetscReal   h, x = 0.0;
360:   PetscScalar uu;

362:   PetscFunctionBegin;
363:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, solution));
364:   h = 1.0 / ((PetscReal)(n + 1));
365:   for (i = 0; i < n; i++) {
366:     x += h;
367:     uu = x * (1. - x);
368:     PetscCall(VecSetValues(*solution, 1, &i, &uu, INSERT_VALUES));
369:   }
370:   PetscFunctionReturn(PETSC_SUCCESS);
371: }

373: PetscErrorCode CalculateError(Vec solution, Vec u, Vec r, PetscReal *e)
374: {
375:   PetscFunctionBegin;
376:   PetscCall(VecNorm(r, NORM_2, e + 2));
377:   PetscCall(VecWAXPY(r, -1.0, u, solution));
378:   PetscCall(VecNorm(r, NORM_2, e));
379:   PetscCall(VecNorm(r, NORM_1, e + 1));
380:   PetscFunctionReturn(PETSC_SUCCESS);
381: }

383: /*TEST

385:    test:

387: TEST*/