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: SLES csles,sles[20],slesmg;
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: KSP ksp,kspmg;
44: PC pcmg,pc;
45: PetscTruth flg;
47: PetscInitialize(&Argc,&Args,(char *)0,help);
49: PetscOptionsGetInt(PETSC_NULL,"-x",&x_mesh,PETSC_NULL);
50: PetscOptionsGetInt(PETSC_NULL,"-l",&levels,PETSC_NULL);
51: PetscOptionsGetInt(PETSC_NULL,"-c",&cycles,PETSC_NULL);
52: PetscOptionsGetInt(PETSC_NULL,"-smooths",&smooths,PETSC_NULL);
53: PetscOptionsHasName(PETSC_NULL,"-a",&flg);
54: if (flg) {am = MGADDITIVE;}
55: PetscOptionsHasName(PETSC_NULL,"-f",&flg);
56: if (flg) {am = MGFULL;}
57: PetscOptionsHasName(PETSC_NULL,"-j",&flg);
58: if (flg) {use_jacobi = 1;}
59:
60: PetscMalloc(levels*sizeof(int),&N);
61: N[0] = x_mesh;
62: for (i=1; i<levels; i++) {
63: N[i] = N[i-1]/2;
64: if (N[i] < 1) {SETERRQ(1,"Too many levels");}
65: }
67: Create1dLaplacian(N[levels-1],&cmat);
69: SLESCreate(PETSC_COMM_WORLD,&slesmg);
70: SLESGetPC(slesmg,&pcmg);
71: SLESGetKSP(slesmg,&kspmg);
72: SLESSetFromOptions(slesmg);
73: PCSetType(pcmg,PCMG);
74: MGSetLevels(pcmg,levels,PETSC_NULL);
75: MGSetType(pcmg,am);
77: MGGetCoarseSolve(pcmg,&csles);
78: SLESSetOperators(csles,cmat,cmat,DIFFERENT_NONZERO_PATTERN);
79: SLESGetPC(csles,&pc);
80: PCSetType(pc,PCLU);
81: SLESGetKSP(csles,&ksp);
82: KSPSetType(ksp,KSPPREONLY);
84: /* zero is finest level */
85: for (i=0; i<levels-1; i++) {
86: MGSetResidual(pcmg,levels - 1 - i,residual,(Mat)0);
87: MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],(void *)0,&mat[i]);
88: MatShellSetOperation(mat[i],MATOP_MULT,(void(*)(void))restrct);
89: MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void(*)(void))interpolate);
90: MGSetInterpolate(pcmg,levels - 1 - i,mat[i]);
91: MGSetRestriction(pcmg,levels - 1 - i,mat[i]);
92: MGSetCyclesOnLevel(pcmg,levels - 1 - i,cycles);
94: /* set smoother */
95: MGGetSmoother(pcmg,levels - 1 - i,&sles[i]);
96: SLESGetPC(sles[i],&pc);
97: PCSetType(pc,PCSHELL);
98: PCShellSetName(pc,"user_precond");
99: PCShellGetName(pc,&shellname);
100: PetscPrintf(PETSC_COMM_WORLD,"level=%d, PCShell name is %s\n",i,shellname);
102: /* this is a dummy! since SLES requires a matrix passed in */
103: SLESSetOperators(sles[i],mat[i],mat[i],DIFFERENT_NONZERO_PATTERN);
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: PCShellSetApplyRichardson(pc,gauss_seidel,(void *)0);
110: if (use_jacobi) {
111: PCShellSetApplyRichardson(pc,jacobi,(void *)0);
112: }
113: SLESGetKSP(sles[i],&ksp);
114: KSPSetType(ksp,KSPRICHARDSON);
115: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
116: KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);
118: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
119: X[levels - 1 - i] = x;
120: MGSetX(pcmg,levels - 1 - i,x);
121: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
122: B[levels -1 - i] = x;
123: MGSetRhs(pcmg,levels - 1 - i,x);
124: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
125: R[levels - 1 - i] = x;
126: MGSetR(pcmg,levels - 1 - i,x);
127: }
128: /* create coarse level vectors */
129: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
130: MGSetX(pcmg,0,x); X[0] = x;
131: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
132: MGSetRhs(pcmg,0,x); B[0] = x;
133: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
134: MGSetR(pcmg,0,x); R[0] = x;
136: /* create matrix multiply for finest level */
137: MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],(void *)0,&fmat);
138: MatShellSetOperation(fmat,MATOP_MULT,(void(*)(void))amult);
139: SLESSetOperators(slesmg,fmat,fmat,DIFFERENT_NONZERO_PATTERN);
141: CalculateSolution(N[0],&solution);
142: CalculateRhs(B[levels-1]);
143: VecSet(&zero,X[levels-1]);
145: if (MGCheck(pcmg)) {SETERRQ(PETSC_ERR_PLIB, "MGCheck failed");}
146:
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,"l_2 error %g max error %g resi %g\n",e[0],e[1],e[2]);
151: SLESSolve(slesmg,B[levels-1],X[levels-1],&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,e[0],e[1],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: VecDestroy(R[i]);
165: }
166: for (i=0; i<levels-1; i++) {
167: MatDestroy(mat[i]);
168: }
169: MatDestroy(cmat);
170: MatDestroy(fmat);
171: SLESDestroy(slesmg);
172: PetscFinalize();
173: return 0;
174: }
176: /* --------------------------------------------------------------------- */
179: int residual(Mat mat,Vec bb,Vec xx,Vec rr)
180: {
181: int i,n1,ierr;
182: PetscScalar *b,*x,*r;
185: VecGetSize(bb,&n1);
186: VecGetArray(bb,&b);
187: VecGetArray(xx,&x);
188: VecGetArray(rr,&r);
189: n1--;
190: r[0] = b[0] + x[1] - 2.0*x[0];
191: r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
192: for (i=1; i<n1; i++) {
193: r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
194: }
195: VecRestoreArray(bb,&b);
196: VecRestoreArray(xx,&x);
197: VecRestoreArray(rr,&r);
198: return(0);
199: }
202: int amult(Mat mat,Vec xx,Vec yy)
203: {
204: int i,n1,ierr;
205: PetscScalar *y,*x;
208: VecGetSize(xx,&n1);
209: VecGetArray(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++) {
215: y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
216: }
217: VecRestoreArray(xx,&x);
218: VecRestoreArray(yy,&y);
219: return(0);
220: }
221: /* --------------------------------------------------------------------- */
224: int gauss_seidel(void *ptr,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal atol,PetscReal dtol,int m)
225: {
226: int i,n1,ierr;
227: PetscScalar *x,*b;
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: int jacobi(void *ptr,Vec bb,Vec xx,Vec w,PetscReal rtol,PetscReal atol,PetscReal dtol,int m)
252: {
253: int i,n,n1,ierr;
254: PetscScalar *r,*b,*x;
257: VecGetSize(bb,&n); n1 = n - 1;
258: VecGetArray(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++) {
265: r[i] = .5*(x[i+1] + x[i-1] + b[i]);
266: }
267: r[n1] = .5*(x[n1-1] + b[n1]);
268: for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
269: }
270: VecRestoreArray(bb,&b);
271: VecRestoreArray(xx,&x);
272: VecRestoreArray(w,&r);
273: return(0);
274: }
275: /*
276: We know for this application that yy and zz are the same
277: */
278: /* --------------------------------------------------------------------- */
281: int interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
282: {
283: int i,n,N,i2,ierr;
284: PetscScalar *x,*y;
287: VecGetSize(yy,&N);
288: VecGetArray(xx,&x);
289: VecGetArray(yy,&y);
290: n = N/2;
291: for (i=0; i<n; i++) {
292: i2 = 2*i;
293: y[i2] += .5*x[i];
294: y[i2+1] += x[i];
295: y[i2+2] += .5*x[i];
296: }
297: VecRestoreArray(xx,&x);
298: VecRestoreArray(yy,&y);
299: return(0);
300: }
301: /* --------------------------------------------------------------------- */
304: int restrct(Mat mat,Vec rr,Vec bb)
305: {
306: int i,n,N,i2,ierr;
307: PetscScalar *r,*b;
310: VecGetSize(rr,&N);
311: VecGetArray(rr,&r);
312: VecGetArray(bb,&b);
313: n = N/2;
315: for (i=0; i<n; i++) {
316: i2 = 2*i;
317: b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
318: }
319: VecRestoreArray(rr,&r);
320: VecRestoreArray(bb,&b);
321: return(0);
322: }
323: /* --------------------------------------------------------------------- */
326: int Create1dLaplacian(int n,Mat *mat)
327: {
328: PetscScalar mone = -1.0,two = 2.0;
329: int ierr,i,idx;
332: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,mat);
333:
334: idx= n-1;
335: MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
336: for (i=0; i<n-1; i++) {
337: MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
338: idx = i+1;
339: MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
340: MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
341: }
342: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
343: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
344: return(0);
345: }
346: /* --------------------------------------------------------------------- */
349: int CalculateRhs(Vec u)
350: {
351: int i,n,ierr;
352: PetscReal h,x = 0.0;
353: PetscScalar uu;
356: VecGetSize(u,&n);
357: h = 1.0/((PetscReal)(n+1));
358: for (i=0; i<n; i++) {
359: x += h; uu = 2.0*h*h;
360: VecSetValues(u,1,&i,&uu,INSERT_VALUES);
361: }
363: return(0);
364: }
365: /* --------------------------------------------------------------------- */
368: int CalculateSolution(int n,Vec *solution)
369: {
370: int i,ierr;
371: PetscReal h,x = 0.0;
372: PetscScalar uu;
375: VecCreateSeq(PETSC_COMM_SELF,n,solution);
376: h = 1.0/((PetscReal)(n+1));
377: for (i=0; i<n; i++) {
378: x += h; uu = x*(1.-x);
379: VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
380: }
381: return(0);
382: }
383: /* --------------------------------------------------------------------- */
386: int CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
387: {
388: PetscScalar mone = -1.0;
389: int ierr;
392: VecNorm(r,NORM_2,e+2);
393: VecWAXPY(&mone,u,solution,r);
394: VecNorm(r,NORM_2,e);
395: VecNorm(r,NORM_1,e+1);
396: return(0);
397: }