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";
10
11 /*
12 This is not a good example to understand the use of multigrid with PETSc.
13 */
14
15 #include <petscksp.h>
16
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);
28
main(int Argc,char ** Args)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;
41
42 PetscCall(PetscInitialize(&Argc, &Args, NULL, 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));
48
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;
54
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 }
61
62 PetscCall(Create1dLaplacian(N[levels - 1], &cmat));
63
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));
70
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));
76
77 /* zero is finest level */
78 for (i = 0; i < levels - 1; i++) {
79 Mat dummy;
80
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, (PetscErrorCodeFn *)restrct));
84 PetscCall(MatShellSetOperation(mat[i], MATOP_MULT_TRANSPOSE_ADD, (PetscErrorCodeFn *)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));
88
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));
96
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, (PetscErrorCodeFn *)amult));
100 PetscCall(KSPSetOperators(ksp[i], dummy, dummy));
101 PetscCall(MatDestroy(&dummy));
102
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_CURRENT, PETSC_CURRENT, PETSC_CURRENT, smooths));
114
115 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x));
116
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));
120
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));
124
125 R[levels - 1 - i] = x;
126
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;
136
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, (PetscErrorCodeFn *)amult));
140 PetscCall(KSPSetOperators(kspmg, fmat, fmat));
141
142 PetscCall(CalculateSolution(N[0], &solution));
143 PetscCall(CalculateRhs(B[levels - 1]));
144 PetscCall(VecSet(X[levels - 1], 0.0));
145
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]));
149
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]));
155
156 PetscCall(PetscFree(N));
157 PetscCall(VecDestroy(&solution));
158
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 }
173
residual(Mat mat,Vec bb,Vec xx,Vec rr)174 PetscErrorCode residual(Mat mat, Vec bb, Vec xx, Vec rr)
175 {
176 PetscInt i, n1;
177 PetscScalar *x, *r;
178 const PetscScalar *b;
179
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 }
194
amult(Mat mat,Vec xx,Vec yy)195 PetscErrorCode amult(Mat mat, Vec xx, Vec yy)
196 {
197 PetscInt i, n1;
198 PetscScalar *y;
199 const PetscScalar *x;
200
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 }
213
apply_pc(PC pc,Vec bb,Vec xx)214 PetscErrorCode apply_pc(PC pc, Vec bb, Vec xx)
215 {
216 PetscFunctionBegin;
217 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented");
218 }
219
gauss_seidel(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscBool guesszero,PetscInt * its,PCRichardsonConvergedReason * reason)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;
225
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 }
244
jacobi_smoother(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscBool guesszero,PetscInt * its,PCRichardsonConvergedReason * reason)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;
250
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));
259
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 */
274
interpolate(Mat mat,Vec xx,Vec yy,Vec zz)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;
280
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 }
296
restrct(Mat mat,Vec rr,Vec bb)297 PetscErrorCode restrct(Mat mat, Vec rr, Vec bb)
298 {
299 PetscInt i, n, N, i2;
300 PetscScalar *b;
301 const PetscScalar *r;
302
303 PetscFunctionBegin;
304 PetscCall(VecGetSize(rr, &N));
305 PetscCall(VecGetArrayRead(rr, &r));
306 PetscCall(VecGetArray(bb, &b));
307 n = N / 2;
308
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 }
317
Create1dLaplacian(PetscInt n,Mat * mat)318 PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat)
319 {
320 PetscFunctionBeginUser;
321 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat));
322 PetscCall(MatSetValue(*mat, n - 1, n - 1, 2.0, INSERT_VALUES));
323 for (PetscInt i = 0; i < n - 1; i++) {
324 PetscCall(MatSetValue(*mat, i, i, 2.0, INSERT_VALUES));
325 PetscCall(MatSetValue(*mat, i + 1, i, -1.0, INSERT_VALUES));
326 PetscCall(MatSetValue(*mat, i, i + 1, -1.0, INSERT_VALUES));
327 }
328 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
329 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
330 PetscFunctionReturn(PETSC_SUCCESS);
331 }
332
CalculateRhs(Vec u)333 PetscErrorCode CalculateRhs(Vec u)
334 {
335 PetscInt i, n;
336 PetscReal h;
337 PetscScalar uu;
338
339 PetscFunctionBegin;
340 PetscCall(VecGetSize(u, &n));
341 h = 1.0 / ((PetscReal)(n + 1));
342 for (i = 0; i < n; i++) {
343 uu = 2.0 * h * h;
344 PetscCall(VecSetValues(u, 1, &i, &uu, INSERT_VALUES));
345 }
346 PetscFunctionReturn(PETSC_SUCCESS);
347 }
348
CalculateSolution(PetscInt n,Vec * solution)349 PetscErrorCode CalculateSolution(PetscInt n, Vec *solution)
350 {
351 PetscInt i;
352 PetscReal h, x = 0.0;
353 PetscScalar uu;
354
355 PetscFunctionBegin;
356 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, solution));
357 h = 1.0 / ((PetscReal)(n + 1));
358 for (i = 0; i < n; i++) {
359 x += h;
360 uu = x * (1. - x);
361 PetscCall(VecSetValues(*solution, 1, &i, &uu, INSERT_VALUES));
362 }
363 PetscFunctionReturn(PETSC_SUCCESS);
364 }
365
CalculateError(Vec solution,Vec u,Vec r,PetscReal * e)366 PetscErrorCode CalculateError(Vec solution, Vec u, Vec r, PetscReal *e)
367 {
368 PetscFunctionBegin;
369 PetscCall(VecNorm(r, NORM_2, e + 2));
370 PetscCall(VecWAXPY(r, -1.0, u, solution));
371 PetscCall(VecNorm(r, NORM_2, e));
372 PetscCall(VecNorm(r, NORM_1, e + 1));
373 PetscFunctionReturn(PETSC_SUCCESS);
374 }
375
376 /*TEST
377
378 test:
379
380 TEST*/
381