Actual source code: ex5.c
petsc-3.6.1 2015-08-06
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: */
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(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);
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;
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;
44: PetscInitialize(&Argc,&Args,(char*)0,help);
46: PetscOptionsGetInt(NULL,"-x",&x_mesh,NULL);
47: PetscOptionsGetInt(NULL,"-l",&levels,NULL);
48: PetscOptionsGetInt(NULL,"-c",&cycles,NULL);
49: PetscOptionsGetInt(NULL,"-smooths",&smooths,NULL);
50: PetscOptionsHasName(NULL,"-a",&flg);
52: if (flg) am = PC_MG_ADDITIVE;
53: PetscOptionsHasName(NULL,"-f",&flg);
54: if (flg) am = PC_MG_FULL;
55: PetscOptionsHasName(NULL,"-j",&flg);
56: if (flg) use_jacobi = 1;
58: PetscMalloc1(levels,&N);
59: N[0] = x_mesh;
60: for (i=1; i<levels; i++) {
61: N[i] = N[i-1]/2;
62: if (N[i] < 1) SETERRQ(PETSC_COMM_WORLD,1,"Too many levels");
63: }
65: Create1dLaplacian(N[levels-1],&cmat);
67: KSPCreate(PETSC_COMM_WORLD,&kspmg);
68: KSPGetPC(kspmg,&pcmg);
69: KSPSetFromOptions(kspmg);
70: PCSetType(pcmg,PCMG);
71: PCMGSetLevels(pcmg,levels,NULL);
72: PCMGSetType(pcmg,am);
74: PCMGGetCoarseSolve(pcmg,&cksp);
75: KSPSetOperators(cksp,cmat,cmat);
76: KSPGetPC(cksp,&pc);
77: PCSetType(pc,PCLU);
78: KSPSetType(cksp,KSPPREONLY);
80: /* zero is finest level */
81: for (i=0; i<levels-1; i++) {
82: PCMGSetResidual(pcmg,levels - 1 - i,residual,(Mat)0);
83: MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],(void*)0,&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: PCMGSetCyclesOnLevel(pcmg,levels - 1 - i,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 a dummy! since KSP requires a matrix passed in */
99: KSPSetOperators(ksp[i],mat[i],mat[i]);
100: /*
101: We override the matrix passed in by forcing it to use Richardson with
102: a user provided application. This is non-standard and this practice
103: should be avoided.
104: */
105: PCShellSetApplyRichardson(pc,gauss_seidel);
106: if (use_jacobi) {
107: PCShellSetApplyRichardson(pc,jacobi);
108: }
109: KSPSetType(ksp[i],KSPRICHARDSON);
110: KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);
111: KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);
113: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
115: X[levels - 1 - i] = x;
116: if (i > 0) {
117: PCMGSetX(pcmg,levels - 1 - i,x);
118: }
119: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
121: B[levels -1 - i] = x;
122: if (i > 0) {
123: PCMGSetRhs(pcmg,levels - 1 - i,x);
124: }
125: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
127: R[levels - 1 - i] = x;
129: PCMGSetR(pcmg,levels - 1 - i,x);
130: }
131: /* create coarse level vectors */
132: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
133: PCMGSetX(pcmg,0,x); X[0] = x;
134: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
135: PCMGSetRhs(pcmg,0,x); B[0] = x;
137: /* create matrix multiply for finest level */
138: MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],(void*)0,&fmat);
139: MatShellSetOperation(fmat,MATOP_MULT,(void (*)(void))amult);
140: KSPSetOperators(kspmg,fmat,fmat);
142: CalculateSolution(N[0],&solution);
143: CalculateRhs(B[levels-1]);
144: VecSet(X[levels-1],0.0);
146: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
147: CalculateError(solution,X[levels-1],R[levels-1],e);
148: PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",(double)e[0],(double)e[1],(double)e[2]);
150: KSPSolve(kspmg,B[levels-1],X[levels-1]);
151: KSPGetIterationNumber(kspmg,&its);
152: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
153: CalculateError(solution,X[levels-1],R[levels-1],e);
154: 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]);
156: PetscFree(N);
157: 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: VecDestroy(&X[i]);
163: VecDestroy(&B[i]);
164: if (i) {VecDestroy(&R[i]);}
165: }
166: for (i=0; i<levels-1; i++) {
167: MatDestroy(&mat[i]);
168: }
169: MatDestroy(&cmat);
170: MatDestroy(&fmat);
171: KSPDestroy(&kspmg);
172: PetscFinalize();
173: return 0;
174: }
176: /* --------------------------------------------------------------------- */
179: PetscErrorCode residual(Mat mat,Vec bb,Vec xx,Vec rr)
180: {
181: PetscInt i,n1;
182: PetscErrorCode ierr;
183: PetscScalar *x,*r;
184: const PetscScalar *b;
187: VecGetSize(bb,&n1);
188: VecGetArrayRead(bb,&b);
189: VecGetArray(xx,&x);
190: VecGetArray(rr,&r);
191: n1--;
192: r[0] = b[0] + x[1] - 2.0*x[0];
193: r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
194: for (i=1; i<n1; i++) r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
195: VecRestoreArrayRead(bb,&b);
196: VecRestoreArray(xx,&x);
197: VecRestoreArray(rr,&r);
198: return(0);
199: }
202: PetscErrorCode amult(Mat mat,Vec xx,Vec yy)
203: {
204: PetscInt i,n1;
205: PetscErrorCode ierr;
206: PetscScalar *y;
207: const PetscScalar *x;
210: VecGetSize(xx,&n1);
211: VecGetArrayRead(xx,&x);
212: VecGetArray(yy,&y);
213: n1--;
214: y[0] = -x[1] + 2.0*x[0];
215: y[n1] = -x[n1-1] + 2.0*x[n1];
216: for (i=1; i<n1; i++) y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
217: VecRestoreArrayRead(xx,&x);
218: VecRestoreArray(yy,&y);
219: return(0);
220: }
221: /* --------------------------------------------------------------------- */
224: 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)
225: {
226: PetscInt i,n1;
227: PetscErrorCode ierr;
228: PetscScalar *x;
229: 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: /* --------------------------------------------------------------------- */
251: PetscErrorCode jacobi(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscBool guesszero,PetscInt *its,PCRichardsonConvergedReason *reason)
252: {
253: PetscInt i,n,n1;
254: PetscErrorCode ierr;
255: PetscScalar *r,*x;
256: const PetscScalar *b;
259: *its = m;
260: *reason = PCRICHARDSON_CONVERGED_ITS;
261: VecGetSize(bb,&n); n1 = n - 1;
262: VecGetArrayRead(bb,&b);
263: VecGetArray(xx,&x);
264: VecGetArray(w,&r);
266: while (m--) {
267: r[0] = .5*(x[1] + b[0]);
268: for (i=1; i<n1; i++) r[i] = .5*(x[i+1] + x[i-1] + b[i]);
269: r[n1] = .5*(x[n1-1] + b[n1]);
270: for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
271: }
272: VecRestoreArrayRead(bb,&b);
273: VecRestoreArray(xx,&x);
274: VecRestoreArray(w,&r);
275: return(0);
276: }
277: /*
278: We know for this application that yy and zz are the same
279: */
280: /* --------------------------------------------------------------------- */
283: PetscErrorCode interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
284: {
285: PetscInt i,n,N,i2;
286: PetscErrorCode ierr;
287: PetscScalar *y;
288: const PetscScalar *x;
291: VecGetSize(yy,&N);
292: VecGetArrayRead(xx,&x);
293: VecGetArray(yy,&y);
294: n = N/2;
295: for (i=0; i<n; i++) {
296: i2 = 2*i;
297: y[i2] += .5*x[i];
298: y[i2+1] += x[i];
299: y[i2+2] += .5*x[i];
300: }
301: VecRestoreArrayRead(xx,&x);
302: VecRestoreArray(yy,&y);
303: return(0);
304: }
305: /* --------------------------------------------------------------------- */
308: PetscErrorCode restrct(Mat mat,Vec rr,Vec bb)
309: {
310: PetscInt i,n,N,i2;
311: PetscErrorCode ierr;
312: PetscScalar *b;
313: const PetscScalar *r;
316: VecGetSize(rr,&N);
317: VecGetArrayRead(rr,&r);
318: VecGetArray(bb,&b);
319: n = N/2;
321: for (i=0; i<n; i++) {
322: i2 = 2*i;
323: b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
324: }
325: VecRestoreArrayRead(rr,&r);
326: VecRestoreArray(bb,&b);
327: return(0);
328: }
329: /* --------------------------------------------------------------------- */
332: PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
333: {
334: PetscScalar mone = -1.0,two = 2.0;
335: PetscInt i,idx;
339: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,mat);
341: idx = n-1;
342: MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
343: for (i=0; i<n-1; i++) {
344: MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
345: idx = i+1;
346: MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
347: MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
348: }
349: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
350: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
351: return(0);
352: }
353: /* --------------------------------------------------------------------- */
356: PetscErrorCode CalculateRhs(Vec u)
357: {
359: PetscInt i,n;
360: PetscReal h,x = 0.0;
361: PetscScalar uu;
364: VecGetSize(u,&n);
365: h = 1.0/((PetscReal)(n+1));
366: for (i=0; i<n; i++) {
367: x += h; uu = 2.0*h*h;
368: VecSetValues(u,1,&i,&uu,INSERT_VALUES);
369: }
370: return(0);
371: }
372: /* --------------------------------------------------------------------- */
375: PetscErrorCode CalculateSolution(PetscInt n,Vec *solution)
376: {
378: PetscInt i;
379: PetscReal h,x = 0.0;
380: PetscScalar uu;
383: VecCreateSeq(PETSC_COMM_SELF,n,solution);
384: h = 1.0/((PetscReal)(n+1));
385: for (i=0; i<n; i++) {
386: x += h; uu = x*(1.-x);
387: VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
388: }
389: return(0);
390: }
391: /* --------------------------------------------------------------------- */
394: PetscErrorCode CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
395: {
399: VecNorm(r,NORM_2,e+2);
400: VecWAXPY(r,-1.0,u,solution);
401: VecNorm(r,NORM_2,e);
402: VecNorm(r,NORM_1,e+1);
403: return(0);
404: }