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