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=%D, 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) {
112: PCShellSetApplyRichardson(pc,jacobi_smoother);
113: }
114: KSPSetType(ksp[i],KSPRICHARDSON);
115: KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);
116: KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);
118: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
120: X[levels - 1 - i] = x;
121: if (i > 0) {
122: PCMGSetX(pcmg,levels - 1 - i,x);
123: }
124: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
126: B[levels -1 - i] = x;
127: if (i > 0) {
128: PCMGSetRhs(pcmg,levels - 1 - i,x);
129: }
130: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
132: R[levels - 1 - i] = x;
134: PCMGSetR(pcmg,levels - 1 - i,x);
135: }
136: /* create coarse level vectors */
137: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
138: PCMGSetX(pcmg,0,x); X[0] = x;
139: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
140: PCMGSetRhs(pcmg,0,x); B[0] = x;
142: /* create matrix multiply for finest level */
143: MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],NULL,&fmat);
144: MatShellSetOperation(fmat,MATOP_MULT,(void (*)(void))amult);
145: KSPSetOperators(kspmg,fmat,fmat);
147: CalculateSolution(N[0],&solution);
148: CalculateRhs(B[levels-1]);
149: VecSet(X[levels-1],0.0);
151: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
152: CalculateError(solution,X[levels-1],R[levels-1],e);
153: PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",(double)e[0],(double)e[1],(double)e[2]);
155: KSPSolve(kspmg,B[levels-1],X[levels-1]);
156: KSPGetIterationNumber(kspmg,&its);
157: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
158: CalculateError(solution,X[levels-1],R[levels-1],e);
159: 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]);
161: PetscFree(N);
162: VecDestroy(&solution);
164: /* note we have to keep a list of all vectors allocated, this is
165: not ideal, but putting it in MGDestroy is not so good either*/
166: for (i=0; i<levels; i++) {
167: VecDestroy(&X[i]);
168: VecDestroy(&B[i]);
169: if (i) VecDestroy(&R[i]);
170: }
171: for (i=0; i<levels-1; i++) {
172: MatDestroy(&mat[i]);
173: }
174: MatDestroy(&cmat);
175: MatDestroy(&fmat);
176: KSPDestroy(&kspmg);
177: PetscFinalize();
178: return 0;
179: }
181: /* --------------------------------------------------------------------- */
182: PetscErrorCode residual(Mat mat,Vec bb,Vec xx,Vec rr)
183: {
184: PetscInt i,n1;
185: PetscScalar *x,*r;
186: const PetscScalar *b;
188: VecGetSize(bb,&n1);
189: VecGetArrayRead(bb,&b);
190: VecGetArray(xx,&x);
191: VecGetArray(rr,&r);
192: n1--;
193: r[0] = b[0] + x[1] - 2.0*x[0];
194: r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
195: for (i=1; i<n1; i++) r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
196: VecRestoreArrayRead(bb,&b);
197: VecRestoreArray(xx,&x);
198: VecRestoreArray(rr,&r);
199: return 0;
200: }
202: PetscErrorCode amult(Mat mat,Vec xx,Vec yy)
203: {
204: PetscInt i,n1;
205: PetscScalar *y;
206: const PetscScalar *x;
208: VecGetSize(xx,&n1);
209: VecGetArrayRead(xx,&x);
210: VecGetArray(yy,&y);
211: n1--;
212: y[0] = -x[1] + 2.0*x[0];
213: y[n1] = -x[n1-1] + 2.0*x[n1];
214: for (i=1; i<n1; i++) y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
215: VecRestoreArrayRead(xx,&x);
216: VecRestoreArray(yy,&y);
217: return 0;
218: }
220: /* --------------------------------------------------------------------- */
221: PetscErrorCode apply_pc(PC pc,Vec bb,Vec xx)
222: {
223: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not implemented");
224: }
226: 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)
227: {
228: PetscInt i,n1;
229: PetscScalar *x;
230: const PetscScalar *b;
232: *its = m;
233: *reason = PCRICHARDSON_CONVERGED_ITS;
234: VecGetSize(bb,&n1); n1--;
235: VecGetArrayRead(bb,&b);
236: VecGetArray(xx,&x);
237: while (m--) {
238: x[0] = .5*(x[1] + b[0]);
239: for (i=1; i<n1; i++) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
240: x[n1] = .5*(x[n1-1] + b[n1]);
241: for (i=n1-1; i>0; i--) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
242: x[0] = .5*(x[1] + b[0]);
243: }
244: VecRestoreArrayRead(bb,&b);
245: VecRestoreArray(xx,&x);
246: return 0;
247: }
248: /* --------------------------------------------------------------------- */
249: 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)
250: {
251: PetscInt i,n,n1;
252: PetscScalar *r,*x;
253: const PetscScalar *b;
255: *its = m;
256: *reason = PCRICHARDSON_CONVERGED_ITS;
257: VecGetSize(bb,&n); n1 = n - 1;
258: VecGetArrayRead(bb,&b);
259: VecGetArray(xx,&x);
260: VecGetArray(w,&r);
262: while (m--) {
263: r[0] = .5*(x[1] + b[0]);
264: for (i=1; i<n1; i++) r[i] = .5*(x[i+1] + x[i-1] + b[i]);
265: r[n1] = .5*(x[n1-1] + b[n1]);
266: for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
267: }
268: VecRestoreArrayRead(bb,&b);
269: VecRestoreArray(xx,&x);
270: VecRestoreArray(w,&r);
271: return 0;
272: }
273: /*
274: We know for this application that yy and zz are the same
275: */
276: /* --------------------------------------------------------------------- */
277: PetscErrorCode interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
278: {
279: PetscInt i,n,N,i2;
280: PetscScalar *y;
281: const PetscScalar *x;
283: VecGetSize(yy,&N);
284: VecGetArrayRead(xx,&x);
285: 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: VecRestoreArrayRead(xx,&x);
294: VecRestoreArray(yy,&y);
295: return 0;
296: }
297: /* --------------------------------------------------------------------- */
298: PetscErrorCode restrct(Mat mat,Vec rr,Vec bb)
299: {
300: PetscInt i,n,N,i2;
301: PetscScalar *b;
302: const PetscScalar *r;
304: VecGetSize(rr,&N);
305: VecGetArrayRead(rr,&r);
306: 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: VecRestoreArrayRead(rr,&r);
314: VecRestoreArray(bb,&b);
315: return 0;
316: }
317: /* --------------------------------------------------------------------- */
318: PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
319: {
320: PetscScalar mone = -1.0,two = 2.0;
321: PetscInt i,idx;
323: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,mat);
325: idx = n-1;
326: MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
327: for (i=0; i<n-1; i++) {
328: MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
329: idx = i+1;
330: MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
331: MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
332: }
333: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
334: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
335: return 0;
336: }
337: /* --------------------------------------------------------------------- */
338: PetscErrorCode CalculateRhs(Vec u)
339: {
340: PetscInt i,n;
341: PetscReal h;
342: PetscScalar uu;
344: VecGetSize(u,&n);
345: h = 1.0/((PetscReal)(n+1));
346: for (i=0; i<n; i++) {
347: uu = 2.0*h*h;
348: VecSetValues(u,1,&i,&uu,INSERT_VALUES);
349: }
350: return 0;
351: }
352: /* --------------------------------------------------------------------- */
353: PetscErrorCode CalculateSolution(PetscInt n,Vec *solution)
354: {
355: PetscInt i;
356: PetscReal h,x = 0.0;
357: PetscScalar uu;
359: VecCreateSeq(PETSC_COMM_SELF,n,solution);
360: h = 1.0/((PetscReal)(n+1));
361: for (i=0; i<n; i++) {
362: x += h; uu = x*(1.-x);
363: VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
364: }
365: return 0;
366: }
367: /* --------------------------------------------------------------------- */
368: PetscErrorCode CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
369: {
370: VecNorm(r,NORM_2,e+2);
371: VecWAXPY(r,-1.0,u,solution);
372: VecNorm(r,NORM_2,e);
373: VecNorm(r,NORM_1,e+1);
374: return 0;
375: }
377: /*TEST
379: test:
381: TEST*/