Actual source code: ex5.c
1: /*$Id: ex5.c,v 1.76 2001/08/07 21:30:37 bsmith Exp $*/
3: static char help[] = "Tests the multigrid code. The input parameters are:\n\
4: -x N Use a mesh in the x direction of N. \n\
5: -c N Use N V-cycles. \n\
6: -l N Use N Levels. \n\
7: -smooths N Use N pre smooths and N post smooths. \n\
8: -j Use Jacobi smoother. \n\
9: -a use additive multigrid \n\
10: -f use full multigrid (preconditioner variant) \n\
11: This example also demonstrates matrix-free methods\n\n";
13: /*
14: This is not a good example to understand the use of multigrid with PETSc.
15: */
16: #include petscmg.h
18: int residual(Mat,Vec,Vec,Vec);
19: int gauss_seidel(void *,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,int);
20: int jacobi(void *,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,int);
21: int interpolate(Mat,Vec,Vec,Vec);
22: int restrct(Mat,Vec,Vec);
23: int Create1dLaplacian(int,Mat*);
24: int CalculateRhs(Vec);
25: int CalculateError(Vec,Vec,Vec,PetscReal*);
26: int CalculateSolution(int,Vec*);
27: int amult(Mat,Vec,Vec);
31: int main(int Argc,char **Args)
32: {
33: int x_mesh = 15,levels = 3,cycles = 1,use_jacobi = 0;
34: int i,smooths = 1,*N;
35: int ierr,its;
36: MGType am = MGMULTIPLICATIVE;
37: Mat cmat,mat[20],fmat;
38: KSP cksp,ksp[20],kspmg;
39: PetscReal e[3]; /* l_2 error,max error, residual */
40: char *shellname;
41: Vec x,solution,X[20],R[20],B[20];
42: PetscScalar zero = 0.0;
43: PC pcmg,pc;
44: PetscTruth flg;
46: PetscInitialize(&Argc,&Args,(char *)0,help);
48: PetscOptionsGetInt(PETSC_NULL,"-x",&x_mesh,PETSC_NULL);
49: PetscOptionsGetInt(PETSC_NULL,"-l",&levels,PETSC_NULL);
50: PetscOptionsGetInt(PETSC_NULL,"-c",&cycles,PETSC_NULL);
51: PetscOptionsGetInt(PETSC_NULL,"-smooths",&smooths,PETSC_NULL);
52: PetscOptionsHasName(PETSC_NULL,"-a",&flg);
53: if (flg) {am = MGADDITIVE;}
54: PetscOptionsHasName(PETSC_NULL,"-f",&flg);
55: if (flg) {am = MGFULL;}
56: PetscOptionsHasName(PETSC_NULL,"-j",&flg);
57: if (flg) {use_jacobi = 1;}
58:
59: PetscMalloc(levels*sizeof(int),&N);
60: N[0] = x_mesh;
61: for (i=1; i<levels; i++) {
62: N[i] = N[i-1]/2;
63: if (N[i] < 1) {SETERRQ(1,"Too many levels");}
64: }
66: Create1dLaplacian(N[levels-1],&cmat);
68: KSPCreate(PETSC_COMM_WORLD,&kspmg);
69: KSPGetPC(kspmg,&pcmg);
70: KSPSetFromOptions(kspmg);
71: PCSetType(pcmg,PCMG);
72: MGSetLevels(pcmg,levels,PETSC_NULL);
73: MGSetType(pcmg,am);
75: MGGetCoarseSolve(pcmg,&cksp);
76: KSPSetOperators(cksp,cmat,cmat,DIFFERENT_NONZERO_PATTERN);
77: KSPGetPC(cksp,&pc);
78: PCSetType(pc,PCLU);
79: KSPSetType(cksp,KSPPREONLY);
81: /* zero is finest level */
82: for (i=0; i<levels-1; i++) {
83: MGSetResidual(pcmg,levels - 1 - i,residual,(Mat)0);
84: MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],(void *)0,&mat[i]);
85: MatShellSetOperation(mat[i],MATOP_MULT,(void(*)(void))restrct);
86: MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void(*)(void))interpolate);
87: MGSetInterpolate(pcmg,levels - 1 - i,mat[i]);
88: MGSetRestriction(pcmg,levels - 1 - i,mat[i]);
89: MGSetCyclesOnLevel(pcmg,levels - 1 - i,cycles);
91: /* set smoother */
92: MGGetSmoother(pcmg,levels - 1 - i,&ksp[i]);
93: KSPGetPC(ksp[i],&pc);
94: PCSetType(pc,PCSHELL);
95: PCShellSetName(pc,"user_precond");
96: PCShellGetName(pc,&shellname);
97: PetscPrintf(PETSC_COMM_WORLD,"level=%d, PCShell name is %s\n",i,shellname);
99: /* this is a dummy! since KSP requires a matrix passed in */
100: KSPSetOperators(ksp[i],mat[i],mat[i],DIFFERENT_NONZERO_PATTERN);
101: /*
102: We override the matrix passed in by forcing it to use Richardson with
103: a user provided application. This is non-standard and this practice
104: should be avoided.
105: */
106: PCShellSetApplyRichardson(pc,gauss_seidel,(void *)0);
107: if (use_jacobi) {
108: PCShellSetApplyRichardson(pc,jacobi,(void *)0);
109: }
110: KSPSetType(ksp[i],KSPRICHARDSON);
111: KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);
112: KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);
114: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
115: X[levels - 1 - i] = x;
116: MGSetX(pcmg,levels - 1 - i,x);
117: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
118: B[levels -1 - i] = x;
119: MGSetRhs(pcmg,levels - 1 - i,x);
120: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
121: R[levels - 1 - i] = x;
122: MGSetR(pcmg,levels - 1 - i,x);
123: }
124: /* create coarse level vectors */
125: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
126: MGSetX(pcmg,0,x); X[0] = x;
127: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
128: MGSetRhs(pcmg,0,x); B[0] = x;
129: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
130: MGSetR(pcmg,0,x); R[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(&zero,X[levels-1]);
141: if (MGCheck(pcmg)) {SETERRQ(PETSC_ERR_PLIB, "MGCheck failed");}
142:
143: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
144: CalculateError(solution,X[levels-1],R[levels-1],e);
145: PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",e[0],e[1],e[2]);
147: KSPSetRhs(kspmg,B[levels-1]);
148: KSPSetSolution(kspmg,X[levels-1]);
149: KSPSolve(kspmg);
150: KSPGetIterationNumber(kspmg,&its);
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,"its %d l_2 error %g max error %g resi %g\n",its,e[0],e[1],e[2]);
155: PetscFree(N);
156: VecDestroy(solution);
158: /* note we have to keep a list of all vectors allocated, this is
159: not ideal, but putting it in MGDestroy is not so good either*/
160: for (i=0; i<levels; i++) {
161: VecDestroy(X[i]);
162: VecDestroy(B[i]);
163: VecDestroy(R[i]);
164: }
165: for (i=0; i<levels-1; i++) {
166: MatDestroy(mat[i]);
167: }
168: MatDestroy(cmat);
169: MatDestroy(fmat);
170: KSPDestroy(kspmg);
171: PetscFinalize();
172: return 0;
173: }
175: /* --------------------------------------------------------------------- */
178: int residual(Mat mat,Vec bb,Vec xx,Vec rr)
179: {
180: int i,n1,ierr;
181: PetscScalar *b,*x,*r;
184: VecGetSize(bb,&n1);
185: VecGetArray(bb,&b);
186: VecGetArray(xx,&x);
187: VecGetArray(rr,&r);
188: n1--;
189: r[0] = b[0] + x[1] - 2.0*x[0];
190: r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
191: for (i=1; i<n1; i++) {
192: r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
193: }
194: VecRestoreArray(bb,&b);
195: VecRestoreArray(xx,&x);
196: VecRestoreArray(rr,&r);
197: return(0);
198: }
201: int amult(Mat mat,Vec xx,Vec yy)
202: {
203: int i,n1,ierr;
204: PetscScalar *y,*x;
207: VecGetSize(xx,&n1);
208: VecGetArray(xx,&x);
209: VecGetArray(yy,&y);
210: n1--;
211: y[0] = -x[1] + 2.0*x[0];
212: y[n1] = -x[n1-1] + 2.0*x[n1];
213: for (i=1; i<n1; i++) {
214: y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
215: }
216: VecRestoreArray(xx,&x);
217: VecRestoreArray(yy,&y);
218: return(0);
219: }
220: /* --------------------------------------------------------------------- */
223: int gauss_seidel(void *ptr,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal atol,PetscReal dtol,int m)
224: {
225: int i,n1,ierr;
226: PetscScalar *x,*b;
229: VecGetSize(bb,&n1); n1--;
230: VecGetArray(bb,&b);
231: VecGetArray(xx,&x);
232: while (m--) {
233: x[0] = .5*(x[1] + b[0]);
234: for (i=1; i<n1; i++) {
235: x[i] = .5*(x[i+1] + x[i-1] + b[i]);
236: }
237: x[n1] = .5*(x[n1-1] + b[n1]);
238: for (i=n1-1; i>0; i--) {
239: x[i] = .5*(x[i+1] + x[i-1] + b[i]);
240: }
241: x[0] = .5*(x[1] + b[0]);
242: }
243: VecRestoreArray(bb,&b);
244: VecRestoreArray(xx,&x);
245: return(0);
246: }
247: /* --------------------------------------------------------------------- */
250: int jacobi(void *ptr,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal atol,PetscReal dtol,int m)
251: {
252: int i,n,n1,ierr;
253: PetscScalar *r,*b,*x;
256: VecGetSize(bb,&n); n1 = n - 1;
257: VecGetArray(bb,&b);
258: VecGetArray(xx,&x);
259: VecGetArray(w,&r);
261: while (m--) {
262: r[0] = .5*(x[1] + b[0]);
263: for (i=1; i<n1; i++) {
264: r[i] = .5*(x[i+1] + x[i-1] + b[i]);
265: }
266: r[n1] = .5*(x[n1-1] + b[n1]);
267: for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
268: }
269: VecRestoreArray(bb,&b);
270: VecRestoreArray(xx,&x);
271: VecRestoreArray(w,&r);
272: return(0);
273: }
274: /*
275: We know for this application that yy and zz are the same
276: */
277: /* --------------------------------------------------------------------- */
280: int interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
281: {
282: int i,n,N,i2,ierr;
283: PetscScalar *x,*y;
286: VecGetSize(yy,&N);
287: VecGetArray(xx,&x);
288: VecGetArray(yy,&y);
289: n = N/2;
290: for (i=0; i<n; i++) {
291: i2 = 2*i;
292: y[i2] += .5*x[i];
293: y[i2+1] += x[i];
294: y[i2+2] += .5*x[i];
295: }
296: VecRestoreArray(xx,&x);
297: VecRestoreArray(yy,&y);
298: return(0);
299: }
300: /* --------------------------------------------------------------------- */
303: int restrct(Mat mat,Vec rr,Vec bb)
304: {
305: int i,n,N,i2,ierr;
306: PetscScalar *r,*b;
309: VecGetSize(rr,&N);
310: VecGetArray(rr,&r);
311: VecGetArray(bb,&b);
312: n = N/2;
314: for (i=0; i<n; i++) {
315: i2 = 2*i;
316: b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
317: }
318: VecRestoreArray(rr,&r);
319: VecRestoreArray(bb,&b);
320: return(0);
321: }
322: /* --------------------------------------------------------------------- */
325: int Create1dLaplacian(int n,Mat *mat)
326: {
327: PetscScalar mone = -1.0,two = 2.0;
328: int ierr,i,idx;
331: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,mat);
332:
333: idx= n-1;
334: MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
335: for (i=0; i<n-1; i++) {
336: MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
337: idx = i+1;
338: MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
339: MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
340: }
341: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
342: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
343: return(0);
344: }
345: /* --------------------------------------------------------------------- */
348: int CalculateRhs(Vec u)
349: {
350: int i,n,ierr;
351: PetscReal h,x = 0.0;
352: PetscScalar uu;
355: VecGetSize(u,&n);
356: h = 1.0/((PetscReal)(n+1));
357: for (i=0; i<n; i++) {
358: x += h; uu = 2.0*h*h;
359: VecSetValues(u,1,&i,&uu,INSERT_VALUES);
360: }
362: return(0);
363: }
364: /* --------------------------------------------------------------------- */
367: int CalculateSolution(int n,Vec *solution)
368: {
369: int i,ierr;
370: PetscReal h,x = 0.0;
371: PetscScalar uu;
374: VecCreateSeq(PETSC_COMM_SELF,n,solution);
375: h = 1.0/((PetscReal)(n+1));
376: for (i=0; i<n; i++) {
377: x += h; uu = x*(1.-x);
378: VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
379: }
380: return(0);
381: }
382: /* --------------------------------------------------------------------- */
385: int CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
386: {
387: PetscScalar mone = -1.0;
388: int ierr;
391: VecNorm(r,NORM_2,e+2);
392: VecWAXPY(&mone,u,solution,r);
393: VecNorm(r,NORM_2,e);
394: VecNorm(r,NORM_1,e+1);
395: return(0);
396: }