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