Actual source code: ex32.c
2: /*
3: Laplacian in 3D. Use for testing BAIJ matrix.
4: Modeled by the partial differential equation
6: - Laplacian u = 1,0 < x,y,z < 1,
8: with boundary conditions
9: u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
10: */
12: static char help[] = "Solves 3D Laplacian using wirebasket based multigrid.\n\n";
14: #include petscda.h
15: #include petscksp.h
16: #include petscmg.h
23: int main(int argc,char **argv)
24: {
26: KSP ksp;
27: PC pc;
28: Vec x,b;
29: DA da;
30: Mat A,Atrans;
31: PetscInt dof=1,M=-8;
32: PetscTruth flg,trans=PETSC_FALSE;
34: PetscInitialize(&argc,&argv,(char *)0,help);
35: PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
36: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
37: PetscOptionsGetTruth(PETSC_NULL,"-trans",&trans,PETSC_NULL);
39: DACreate(PETSC_COMM_WORLD,&da);
40: DASetDim(da,3);
41: DASetPeriodicity(da,DA_NONPERIODIC);
42: DASetStencilType(da,DA_STENCIL_STAR);
43: DASetSizes(da,M,M,M);
44: DASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
45: DASetDof(da,dof);
46: DASetStencilWidth(da,1);
47: DASetVertexDivision(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);
48: DASetFromOptions(da);
50: DACreateGlobalVector(da,&x);
51: DACreateGlobalVector(da,&b);
52: ComputeRHS(da,b);
53: DAGetMatrix(da,MATBAIJ,&A);
54: ComputeMatrix(da,A);
57: /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
58: MatTranspose(A,MAT_INITIAL_MATRIX,&Atrans);
59: MatAXPY(A,1.0,Atrans,DIFFERENT_NONZERO_PATTERN);
60: MatScale(A,0.5);
61: MatDestroy(Atrans);
63: /* Test sbaij matrix */
64: flg = PETSC_FALSE;
65: PetscOptionsGetTruth(PETSC_NULL, "-test_sbaij1", &flg,PETSC_NULL);
66: if (flg){
67: Mat sA;
68: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
69: MatDestroy(A);
70: A = sA;
71: }
73: KSPCreate(PETSC_COMM_WORLD,&ksp);
74: KSPSetFromOptions(ksp);
75: KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
76: KSPGetPC(ksp,&pc);
77: PCSetDA(pc,da);
78:
79: if (trans) {
80: KSPSolveTranspose(ksp,b,x);
81: } else {
82: KSPSolve(ksp,b,x);
83: }
85: /* check final residual */
86: flg = PETSC_FALSE;
87: PetscOptionsGetTruth(PETSC_NULL, "-check_final_residual", &flg,PETSC_NULL);
88: if (flg){
89: Vec b1;
90: PetscReal norm;
91: KSPGetSolution(ksp,&x);
92: VecDuplicate(b,&b1);
93: MatMult(A,x,b1);
94: VecAXPY(b1,-1.0,b);
95: VecNorm(b1,NORM_2,&norm);
96: PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
97: VecDestroy(b1);
98: }
99:
100: KSPDestroy(ksp);
101: VecDestroy(x);
102: VecDestroy(b);
103: MatDestroy(A);
104: DADestroy(da);
105: PetscFinalize();
106: return 0;
107: }
111: PetscErrorCode ComputeRHS(DA da,Vec b)
112: {
114: PetscInt mx,my,mz;
115: PetscScalar h;
118: DAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0);
119: h = 1.0/((mx-1)*(my-1)*(mz-1));
120: VecSet(b,h);
121: return(0);
122: }
123:
126: PetscErrorCode ComputeMatrix(DA da,Mat B)
127: {
129: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
130: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
131: MatStencil row,col;
132:
134: DAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);
135: /* For simplicity, this example only works on mx=my=mz */
136: if ( mx != my || mx != mz) SETERRQ3(1,"This example only works with mx %d = my %d = mz %d\n",mx,my,mz);
138: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
139: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
141: PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
142: v_neighbor = v + dof*dof;
143: PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
144: k3 = 0;
145: for (k1=0; k1<dof; k1++){
146: for (k2=0; k2<dof; k2++){
147: if (k1 == k2){
148: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
149: v_neighbor[k3] = -HxHydHz;
150: } else {
151: v[k3] = k1/(dof*dof); ;
152: v_neighbor[k3] = k2/(dof*dof);
153: }
154: k3++;
155: }
156: }
157: DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
158:
159: for (k=zs; k<zs+zm; k++){
160: for (j=ys; j<ys+ym; j++){
161: for(i=xs; i<xs+xm; i++){
162: row.i = i; row.j = j; row.k = k;
163: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){ /* boudary points */
164: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
165: } else { /* interior points */
166: /* center */
167: col.i = i; col.j = j; col.k = k;
168: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
169:
170: /* x neighbors */
171: col.i = i-1; col.j = j; col.k = k;
172: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
173: col.i = i+1; col.j = j; col.k = k;
174: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
175:
176: /* y neighbors */
177: col.i = i; col.j = j-1; col.k = k;
178: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
179: col.i = i; col.j = j+1; col.k = k;
180: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
181:
182: /* z neighbors */
183: col.i = i; col.j = j; col.k = k-1;
184: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
185: col.i = i; col.j = j; col.k = k+1;
186: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
187: }
188: }
189: }
190: }
191: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
192: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
193: PetscFree(v);
194: return(0);
195: }