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*/