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: */
15: #include petscmg.h
17: PetscErrorCode residual(Mat,Vec,Vec,Vec);
18: PetscErrorCode gauss_seidel(PC,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,PetscInt,PetscTruth,PetscInt*,PCRichardsonConvergedReason*);
19: PetscErrorCode jacobi(PC,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,PetscInt,PetscTruth,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;
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: char *shellname;
40: Vec x,solution,X[20],R[20],B[20];
41: PC pcmg,pc;
42: PetscTruth flg;
44: PetscInitialize(&Argc,&Args,(char *)0,help);
46: PetscOptionsGetInt(PETSC_NULL,"-x",&x_mesh,PETSC_NULL);
47: PetscOptionsGetInt(PETSC_NULL,"-l",&levels,PETSC_NULL);
48: PetscOptionsGetInt(PETSC_NULL,"-c",&cycles,PETSC_NULL);
49: PetscOptionsGetInt(PETSC_NULL,"-smooths",&smooths,PETSC_NULL);
50: PetscOptionsHasName(PETSC_NULL,"-a",&flg);
51: if (flg) {am = PC_MG_ADDITIVE;}
52: PetscOptionsHasName(PETSC_NULL,"-f",&flg);
53: if (flg) {am = PC_MG_FULL;}
54: PetscOptionsHasName(PETSC_NULL,"-j",&flg);
55: if (flg) {use_jacobi = 1;}
56:
57: PetscMalloc(levels*sizeof(PetscInt),&N);
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(1,"Too many levels");}
62: }
64: Create1dLaplacian(N[levels-1],&cmat);
66: KSPCreate(PETSC_COMM_WORLD,&kspmg);
67: KSPGetPC(kspmg,&pcmg);
68: KSPSetFromOptions(kspmg);
69: PCSetType(pcmg,PCMG);
70: PCMGSetLevels(pcmg,levels,PETSC_NULL);
71: PCMGSetType(pcmg,am);
73: PCMGGetCoarseSolve(pcmg,&cksp);
74: KSPSetOperators(cksp,cmat,cmat,DIFFERENT_NONZERO_PATTERN);
75: KSPGetPC(cksp,&pc);
76: PCSetType(pc,PCLU);
77: KSPSetType(cksp,KSPPREONLY);
79: /* zero is finest level */
80: for (i=0; i<levels-1; i++) {
81: PCMGSetResidual(pcmg,levels - 1 - i,residual,(Mat)0);
82: MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],(void*)0,&mat[i]);
83: MatShellSetOperation(mat[i],MATOP_MULT,(void(*)(void))restrct);
84: MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void(*)(void))interpolate);
85: PCMGSetInterpolation(pcmg,levels - 1 - i,mat[i]);
86: PCMGSetRestriction(pcmg,levels - 1 - i,mat[i]);
87: PCMGSetCyclesOnLevel(pcmg,levels - 1 - i,cycles);
89: /* set smoother */
90: PCMGGetSmoother(pcmg,levels - 1 - i,&ksp[i]);
91: KSPGetPC(ksp[i],&pc);
92: PCSetType(pc,PCSHELL);
93: PCShellSetName(pc,"user_precond");
94: PCShellGetName(pc,&shellname);
95: PetscPrintf(PETSC_COMM_WORLD,"level=%D, PCShell name is %s\n",i,shellname);
97: /* this is a dummy! since KSP requires a matrix passed in */
98: KSPSetOperators(ksp[i],mat[i],mat[i],DIFFERENT_NONZERO_PATTERN);
99: /*
100: We override the matrix passed in by forcing it to use Richardson with
101: a user provided application. This is non-standard and this practice
102: should be avoided.
103: */
104: PCShellSetApplyRichardson(pc,gauss_seidel);
105: if (use_jacobi) {
106: PCShellSetApplyRichardson(pc,jacobi);
107: }
108: KSPSetType(ksp[i],KSPRICHARDSON);
109: KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);
110: KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);
112: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
113: X[levels - 1 - i] = x;
114: if (i > 0) {
115: PCMGSetX(pcmg,levels - 1 - i,x);
116: }
117: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
118: B[levels -1 - i] = x;
119: if (i > 0) {
120: PCMGSetRhs(pcmg,levels - 1 - i,x);
121: }
122: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
123: R[levels - 1 - i] = x;
124: PCMGSetR(pcmg,levels - 1 - i,x);
125: }
126: /* create coarse level vectors */
127: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
128: PCMGSetX(pcmg,0,x); X[0] = x;
129: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
130: PCMGSetRhs(pcmg,0,x); B[0] = x;
132: /* create matrix multiply for finest level */
133: MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],(void*)0,&fmat);
134: MatShellSetOperation(fmat,MATOP_MULT,(void(*)(void))amult);
135: KSPSetOperators(kspmg,fmat,fmat,DIFFERENT_NONZERO_PATTERN);
137: CalculateSolution(N[0],&solution);
138: CalculateRhs(B[levels-1]);
139: VecSet(X[levels-1],0.0);
141: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
142: CalculateError(solution,X[levels-1],R[levels-1],e);
143: PetscPrintf(PETSC_COMM_SELF,"l_2 error %G max error %G resi %G\n",e[0],e[1],e[2]);
145: KSPSolve(kspmg,B[levels-1],X[levels-1]);
146: KSPGetIterationNumber(kspmg,&its);
147: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
148: CalculateError(solution,X[levels-1],R[levels-1],e);
149: PetscPrintf(PETSC_COMM_SELF,"its %D l_2 error %G max error %G resi %G\n",its,e[0],e[1],e[2]);
151: PetscFree(N);
152: VecDestroy(solution);
154: /* note we have to keep a list of all vectors allocated, this is
155: not ideal, but putting it in MGDestroy is not so good either*/
156: for (i=0; i<levels; i++) {
157: VecDestroy(X[i]);
158: VecDestroy(B[i]);
159: if(i){VecDestroy(R[i]);}
160: }
161: for (i=0; i<levels-1; i++) {
162: MatDestroy(mat[i]);
163: }
164: MatDestroy(cmat);
165: MatDestroy(fmat);
166: KSPDestroy(kspmg);
167: PetscFinalize();
168: return 0;
169: }
171: /* --------------------------------------------------------------------- */
174: PetscErrorCode residual(Mat mat,Vec bb,Vec xx,Vec rr)
175: {
176: PetscInt i,n1;
178: PetscScalar *b,*x,*r;
181: VecGetSize(bb,&n1);
182: VecGetArray(bb,&b);
183: VecGetArray(xx,&x);
184: 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++) {
189: r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
190: }
191: VecRestoreArray(bb,&b);
192: VecRestoreArray(xx,&x);
193: VecRestoreArray(rr,&r);
194: return(0);
195: }
198: PetscErrorCode amult(Mat mat,Vec xx,Vec yy)
199: {
200: PetscInt i,n1;
202: PetscScalar *y,*x;
205: VecGetSize(xx,&n1);
206: VecGetArray(xx,&x);
207: VecGetArray(yy,&y);
208: n1--;
209: y[0] = -x[1] + 2.0*x[0];
210: y[n1] = -x[n1-1] + 2.0*x[n1];
211: for (i=1; i<n1; i++) {
212: y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
213: }
214: VecRestoreArray(xx,&x);
215: VecRestoreArray(yy,&y);
216: return(0);
217: }
218: /* --------------------------------------------------------------------- */
221: PetscErrorCode gauss_seidel(PC pc,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt m,PetscTruth guesszero,PetscInt *its,PCRichardsonConvergedReason *reason)
222: {
223: PetscInt i,n1;
225: PetscScalar *x,*b;
228: *its = m;
229: *reason = PCRICHARDSON_CONVERGED_ITS;
230: VecGetSize(bb,&n1); n1--;
231: VecGetArray(bb,&b);
232: VecGetArray(xx,&x);
233: while (m--) {
234: x[0] = .5*(x[1] + b[0]);
235: for (i=1; i<n1; i++) {
236: x[i] = .5*(x[i+1] + x[i-1] + b[i]);
237: }
238: x[n1] = .5*(x[n1-1] + b[n1]);
239: for (i=n1-1; i>0; i--) {
240: x[i] = .5*(x[i+1] + x[i-1] + b[i]);
241: }
242: x[0] = .5*(x[1] + b[0]);
243: }
244: VecRestoreArray(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,PetscTruth guesszero,PetscInt *its,PCRichardsonConvergedReason *reason)
252: {
253: PetscInt i,n,n1;
255: PetscScalar *r,*b,*x;
258: *its = m;
259: *reason = PCRICHARDSON_CONVERGED_ITS;
260: VecGetSize(bb,&n); n1 = n - 1;
261: VecGetArray(bb,&b);
262: VecGetArray(xx,&x);
263: VecGetArray(w,&r);
265: while (m--) {
266: r[0] = .5*(x[1] + b[0]);
267: for (i=1; i<n1; i++) {
268: r[i] = .5*(x[i+1] + x[i-1] + b[i]);
269: }
270: r[n1] = .5*(x[n1-1] + b[n1]);
271: for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
272: }
273: VecRestoreArray(bb,&b);
274: VecRestoreArray(xx,&x);
275: VecRestoreArray(w,&r);
276: return(0);
277: }
278: /*
279: We know for this application that yy and zz are the same
280: */
281: /* --------------------------------------------------------------------- */
284: PetscErrorCode interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
285: {
286: PetscInt i,n,N,i2;
288: PetscScalar *x,*y;
291: VecGetSize(yy,&N);
292: VecGetArray(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: VecRestoreArray(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;
312: PetscScalar *r,*b;
315: VecGetSize(rr,&N);
316: VecGetArray(rr,&r);
317: VecGetArray(bb,&b);
318: n = N/2;
320: for (i=0; i<n; i++) {
321: i2 = 2*i;
322: b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
323: }
324: VecRestoreArray(rr,&r);
325: VecRestoreArray(bb,&b);
326: return(0);
327: }
328: /* --------------------------------------------------------------------- */
331: PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
332: {
333: PetscScalar mone = -1.0,two = 2.0;
334: PetscInt i,idx;
338: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,mat);
339:
340: idx= n-1;
341: MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
342: for (i=0; i<n-1; i++) {
343: MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
344: idx = i+1;
345: MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
346: MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
347: }
348: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
349: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
350: return(0);
351: }
352: /* --------------------------------------------------------------------- */
355: PetscErrorCode CalculateRhs(Vec u)
356: {
358: PetscInt i,n;
359: PetscReal h,x = 0.0;
360: PetscScalar uu;
363: VecGetSize(u,&n);
364: h = 1.0/((PetscReal)(n+1));
365: for (i=0; i<n; i++) {
366: x += h; uu = 2.0*h*h;
367: VecSetValues(u,1,&i,&uu,INSERT_VALUES);
368: }
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: }