Actual source code: ex5.c
1: static char help[] = "Tests the multigrid code. The input parameters are:\n\
2: -x N Use a mesh in the x direction of N. \n\
3: -c N Use N V-cycles. \n\
4: -l N Use N Levels. \n\
5: -smooths N Use N pre smooths and N post smooths. \n\
6: -j Use Jacobi smoother. \n\
7: -a use additive multigrid \n\
8: -f use full multigrid (preconditioner variant) \n\
9: This example also demonstrates matrix-free methods\n\n";
11: /*
12: This is not a good example to understand the use of multigrid with PETSc.
13: */
15: #include <petscksp.h>
17: PetscErrorCode residual(Mat, Vec, Vec, Vec);
18: PetscErrorCode gauss_seidel(PC, Vec, Vec, Vec, PetscReal, PetscReal, PetscReal, PetscInt, PetscBool, PetscInt *, PCRichardsonConvergedReason *);
19: PetscErrorCode jacobi_smoother(PC, Vec, Vec, Vec, PetscReal, PetscReal, PetscReal, PetscInt, PetscBool, 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);
27: PetscErrorCode apply_pc(PC, Vec, Vec);
29: int main(int Argc, char **Args)
30: {
31: PetscInt x_mesh = 15, levels = 3, cycles = 1, use_jacobi = 0;
32: PetscInt i, smooths = 1, *N, its;
33: PCMGType am = PC_MG_MULTIPLICATIVE;
34: Mat cmat, mat[20], fmat;
35: KSP cksp, ksp[20], kspmg;
36: PetscReal e[3]; /* l_2 error,max error, residual */
37: const char *shellname;
38: Vec x, solution, X[20], R[20], B[20];
39: PC pcmg, pc;
40: PetscBool flg;
42: PetscCall(PetscInitialize(&Argc, &Args, (char *)0, help));
43: PetscCall(PetscOptionsGetInt(NULL, NULL, "-x", &x_mesh, NULL));
44: PetscCall(PetscOptionsGetInt(NULL, NULL, "-l", &levels, NULL));
45: PetscCall(PetscOptionsGetInt(NULL, NULL, "-c", &cycles, NULL));
46: PetscCall(PetscOptionsGetInt(NULL, NULL, "-smooths", &smooths, NULL));
47: PetscCall(PetscOptionsHasName(NULL, NULL, "-a", &flg));
49: if (flg) am = PC_MG_ADDITIVE;
50: PetscCall(PetscOptionsHasName(NULL, NULL, "-f", &flg));
51: if (flg) am = PC_MG_FULL;
52: PetscCall(PetscOptionsHasName(NULL, NULL, "-j", &flg));
53: if (flg) use_jacobi = 1;
55: PetscCall(PetscMalloc1(levels, &N));
56: N[0] = x_mesh;
57: for (i = 1; i < levels; i++) {
58: N[i] = N[i - 1] / 2;
59: PetscCheck(N[i] >= 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "Too many levels or N is not large enough");
60: }
62: PetscCall(Create1dLaplacian(N[levels - 1], &cmat));
64: PetscCall(KSPCreate(PETSC_COMM_WORLD, &kspmg));
65: PetscCall(KSPGetPC(kspmg, &pcmg));
66: PetscCall(KSPSetFromOptions(kspmg));
67: PetscCall(PCSetType(pcmg, PCMG));
68: PetscCall(PCMGSetLevels(pcmg, levels, NULL));
69: PetscCall(PCMGSetType(pcmg, am));
71: PetscCall(PCMGGetCoarseSolve(pcmg, &cksp));
72: PetscCall(KSPSetOperators(cksp, cmat, cmat));
73: PetscCall(KSPGetPC(cksp, &pc));
74: PetscCall(PCSetType(pc, PCLU));
75: PetscCall(KSPSetType(cksp, KSPPREONLY));
77: /* zero is finest level */
78: for (i = 0; i < levels - 1; i++) {
79: Mat dummy;
81: PetscCall(PCMGSetResidual(pcmg, levels - 1 - i, residual, NULL));
82: PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[i + 1], N[i], N[i + 1], N[i], NULL, &mat[i]));
83: PetscCall(MatShellSetOperation(mat[i], MATOP_MULT, (void (*)(void))restrct));
84: PetscCall(MatShellSetOperation(mat[i], MATOP_MULT_TRANSPOSE_ADD, (void (*)(void))interpolate));
85: PetscCall(PCMGSetInterpolation(pcmg, levels - 1 - i, mat[i]));
86: PetscCall(PCMGSetRestriction(pcmg, levels - 1 - i, mat[i]));
87: PetscCall(PCMGSetCycleTypeOnLevel(pcmg, levels - 1 - i, (PCMGCycleType)cycles));
89: /* set smoother */
90: PetscCall(PCMGGetSmoother(pcmg, levels - 1 - i, &ksp[i]));
91: PetscCall(KSPGetPC(ksp[i], &pc));
92: PetscCall(PCSetType(pc, PCSHELL));
93: PetscCall(PCShellSetName(pc, "user_precond"));
94: PetscCall(PCShellGetName(pc, &shellname));
95: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "level=%" PetscInt_FMT ", PCShell name is %s\n", i, shellname));
97: /* this is not used unless different options are passed to the solver */
98: PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[i], N[i], N[i], N[i], NULL, &dummy));
99: PetscCall(MatShellSetOperation(dummy, MATOP_MULT, (void (*)(void))amult));
100: PetscCall(KSPSetOperators(ksp[i], dummy, dummy));
101: PetscCall(MatDestroy(&dummy));
103: /*
104: We override the matrix passed in by forcing it to use Richardson with
105: a user provided application. This is non-standard and this practice
106: should be avoided.
107: */
108: PetscCall(PCShellSetApply(pc, apply_pc));
109: PetscCall(PCShellSetApplyRichardson(pc, gauss_seidel));
110: if (use_jacobi) PetscCall(PCShellSetApplyRichardson(pc, jacobi_smoother));
111: PetscCall(KSPSetType(ksp[i], KSPRICHARDSON));
112: PetscCall(KSPSetInitialGuessNonzero(ksp[i], PETSC_TRUE));
113: PetscCall(KSPSetTolerances(ksp[i], PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, smooths));
115: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x));
117: X[levels - 1 - i] = x;
118: if (i > 0) PetscCall(PCMGSetX(pcmg, levels - 1 - i, x));
119: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x));
121: B[levels - 1 - i] = x;
122: if (i > 0) PetscCall(PCMGSetRhs(pcmg, levels - 1 - i, x));
123: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x));
125: R[levels - 1 - i] = x;
127: PetscCall(PCMGSetR(pcmg, levels - 1 - i, x));
128: }
129: /* create coarse level vectors */
130: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[levels - 1], &x));
131: PetscCall(PCMGSetX(pcmg, 0, x));
132: X[0] = x;
133: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[levels - 1], &x));
134: PetscCall(PCMGSetRhs(pcmg, 0, x));
135: B[0] = x;
137: /* create matrix multiply for finest level */
138: PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[0], N[0], N[0], N[0], NULL, &fmat));
139: PetscCall(MatShellSetOperation(fmat, MATOP_MULT, (void (*)(void))amult));
140: PetscCall(KSPSetOperators(kspmg, fmat, fmat));
142: PetscCall(CalculateSolution(N[0], &solution));
143: PetscCall(CalculateRhs(B[levels - 1]));
144: PetscCall(VecSet(X[levels - 1], 0.0));
146: PetscCall(residual((Mat)0, B[levels - 1], X[levels - 1], R[levels - 1]));
147: PetscCall(CalculateError(solution, X[levels - 1], R[levels - 1], e));
148: PetscCall(PetscPrintf(PETSC_COMM_SELF, "l_2 error %g max error %g resi %g\n", (double)e[0], (double)e[1], (double)e[2]));
150: PetscCall(KSPSolve(kspmg, B[levels - 1], X[levels - 1]));
151: PetscCall(KSPGetIterationNumber(kspmg, &its));
152: PetscCall(residual((Mat)0, B[levels - 1], X[levels - 1], R[levels - 1]));
153: PetscCall(CalculateError(solution, X[levels - 1], R[levels - 1], e));
154: 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]));
156: PetscCall(PetscFree(N));
157: PetscCall(VecDestroy(&solution));
159: /* note we have to keep a list of all vectors allocated, this is
160: not ideal, but putting it in MGDestroy is not so good either*/
161: for (i = 0; i < levels; i++) {
162: PetscCall(VecDestroy(&X[i]));
163: PetscCall(VecDestroy(&B[i]));
164: if (i) PetscCall(VecDestroy(&R[i]));
165: }
166: for (i = 0; i < levels - 1; i++) PetscCall(MatDestroy(&mat[i]));
167: PetscCall(MatDestroy(&cmat));
168: PetscCall(MatDestroy(&fmat));
169: PetscCall(KSPDestroy(&kspmg));
170: PetscCall(PetscFinalize());
171: return 0;
172: }
174: PetscErrorCode residual(Mat mat, Vec bb, Vec xx, Vec rr)
175: {
176: PetscInt i, n1;
177: PetscScalar *x, *r;
178: const PetscScalar *b;
180: PetscFunctionBegin;
181: PetscCall(VecGetSize(bb, &n1));
182: PetscCall(VecGetArrayRead(bb, &b));
183: PetscCall(VecGetArray(xx, &x));
184: PetscCall(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: PetscCall(VecRestoreArrayRead(bb, &b));
190: PetscCall(VecRestoreArray(xx, &x));
191: PetscCall(VecRestoreArray(rr, &r));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: PetscErrorCode amult(Mat mat, Vec xx, Vec yy)
196: {
197: PetscInt i, n1;
198: PetscScalar *y;
199: const PetscScalar *x;
201: PetscFunctionBegin;
202: PetscCall(VecGetSize(xx, &n1));
203: PetscCall(VecGetArrayRead(xx, &x));
204: PetscCall(VecGetArray(yy, &y));
205: n1--;
206: y[0] = -x[1] + 2.0 * x[0];
207: y[n1] = -x[n1 - 1] + 2.0 * x[n1];
208: for (i = 1; i < n1; i++) y[i] = -x[i + 1] - x[i - 1] + 2.0 * x[i];
209: PetscCall(VecRestoreArrayRead(xx, &x));
210: PetscCall(VecRestoreArray(yy, &y));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: PetscErrorCode apply_pc(PC pc, Vec bb, Vec xx)
215: {
216: PetscFunctionBegin;
217: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented");
218: }
220: 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)
221: {
222: PetscInt i, n1;
223: PetscScalar *x;
224: const PetscScalar *b;
226: PetscFunctionBegin;
227: *its = m;
228: *reason = PCRICHARDSON_CONVERGED_ITS;
229: PetscCall(VecGetSize(bb, &n1));
230: n1--;
231: PetscCall(VecGetArrayRead(bb, &b));
232: PetscCall(VecGetArray(xx, &x));
233: while (m--) {
234: x[0] = .5 * (x[1] + b[0]);
235: for (i = 1; i < n1; i++) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
236: x[n1] = .5 * (x[n1 - 1] + b[n1]);
237: for (i = n1 - 1; i > 0; i--) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
238: x[0] = .5 * (x[1] + b[0]);
239: }
240: PetscCall(VecRestoreArrayRead(bb, &b));
241: PetscCall(VecRestoreArray(xx, &x));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: 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)
246: {
247: PetscInt i, n, n1;
248: PetscScalar *r, *x;
249: const PetscScalar *b;
251: PetscFunctionBegin;
252: *its = m;
253: *reason = PCRICHARDSON_CONVERGED_ITS;
254: PetscCall(VecGetSize(bb, &n));
255: n1 = n - 1;
256: PetscCall(VecGetArrayRead(bb, &b));
257: PetscCall(VecGetArray(xx, &x));
258: PetscCall(VecGetArray(w, &r));
260: while (m--) {
261: r[0] = .5 * (x[1] + b[0]);
262: for (i = 1; i < n1; i++) r[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
263: r[n1] = .5 * (x[n1 - 1] + b[n1]);
264: for (i = 0; i < n; i++) x[i] = (2.0 * r[i] + x[i]) / 3.0;
265: }
266: PetscCall(VecRestoreArrayRead(bb, &b));
267: PetscCall(VecRestoreArray(xx, &x));
268: PetscCall(VecRestoreArray(w, &r));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
271: /*
272: We know for this application that yy and zz are the same
273: */
275: PetscErrorCode interpolate(Mat mat, Vec xx, Vec yy, Vec zz)
276: {
277: PetscInt i, n, N, i2;
278: PetscScalar *y;
279: const PetscScalar *x;
281: PetscFunctionBegin;
282: PetscCall(VecGetSize(yy, &N));
283: PetscCall(VecGetArrayRead(xx, &x));
284: PetscCall(VecGetArray(yy, &y));
285: n = N / 2;
286: for (i = 0; i < n; i++) {
287: i2 = 2 * i;
288: y[i2] += .5 * x[i];
289: y[i2 + 1] += x[i];
290: y[i2 + 2] += .5 * x[i];
291: }
292: PetscCall(VecRestoreArrayRead(xx, &x));
293: PetscCall(VecRestoreArray(yy, &y));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: PetscErrorCode restrct(Mat mat, Vec rr, Vec bb)
298: {
299: PetscInt i, n, N, i2;
300: PetscScalar *b;
301: const PetscScalar *r;
303: PetscFunctionBegin;
304: PetscCall(VecGetSize(rr, &N));
305: PetscCall(VecGetArrayRead(rr, &r));
306: PetscCall(VecGetArray(bb, &b));
307: n = N / 2;
309: for (i = 0; i < n; i++) {
310: i2 = 2 * i;
311: b[i] = (r[i2] + 2.0 * r[i2 + 1] + r[i2 + 2]);
312: }
313: PetscCall(VecRestoreArrayRead(rr, &r));
314: PetscCall(VecRestoreArray(bb, &b));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat)
319: {
320: PetscScalar mone = -1.0, two = 2.0;
321: PetscInt i, idx;
323: PetscFunctionBegin;
324: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat));
326: idx = n - 1;
327: PetscCall(MatSetValues(*mat, 1, &idx, 1, &idx, &two, INSERT_VALUES));
328: for (i = 0; i < n - 1; i++) {
329: PetscCall(MatSetValues(*mat, 1, &i, 1, &i, &two, INSERT_VALUES));
330: idx = i + 1;
331: PetscCall(MatSetValues(*mat, 1, &idx, 1, &i, &mone, INSERT_VALUES));
332: PetscCall(MatSetValues(*mat, 1, &i, 1, &idx, &mone, INSERT_VALUES));
333: }
334: PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
335: PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: PetscErrorCode CalculateRhs(Vec u)
340: {
341: PetscInt i, n;
342: PetscReal h;
343: PetscScalar uu;
345: PetscFunctionBegin;
346: PetscCall(VecGetSize(u, &n));
347: h = 1.0 / ((PetscReal)(n + 1));
348: for (i = 0; i < n; i++) {
349: uu = 2.0 * h * h;
350: PetscCall(VecSetValues(u, 1, &i, &uu, INSERT_VALUES));
351: }
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: PetscErrorCode CalculateSolution(PetscInt n, Vec *solution)
356: {
357: PetscInt i;
358: PetscReal h, x = 0.0;
359: PetscScalar uu;
361: PetscFunctionBegin;
362: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, solution));
363: h = 1.0 / ((PetscReal)(n + 1));
364: for (i = 0; i < n; i++) {
365: x += h;
366: uu = x * (1. - x);
367: PetscCall(VecSetValues(*solution, 1, &i, &uu, INSERT_VALUES));
368: }
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: PetscErrorCode CalculateError(Vec solution, Vec u, Vec r, PetscReal *e)
373: {
374: PetscFunctionBegin;
375: PetscCall(VecNorm(r, NORM_2, e + 2));
376: PetscCall(VecWAXPY(r, -1.0, u, solution));
377: PetscCall(VecNorm(r, NORM_2, e));
378: PetscCall(VecNorm(r, NORM_1, e + 1));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /*TEST
384: test:
386: TEST*/