Actual source code: ex32.c
petsc-3.3-p6 2013-02-11
1: /*
2: Laplacian in 3D. Use for testing BAIJ matrix.
3: Modeled by the partial differential equation
5: - Laplacian u = 1,0 < x,y,z < 1,
7: with boundary conditions
8: u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
9: */
11: static char help[] = "Solves 3D Laplacian using wirebasket based multigrid.\n\n";
13: #include <petscdmda.h>
14: #include <petscksp.h>
15: #include <petscpcmg.h>
17: extern PetscErrorCode ComputeMatrix(DM,Mat);
18: extern PetscErrorCode ComputeRHS(DM,Vec);
22: int main(int argc,char **argv)
23: {
25: KSP ksp;
26: PC pc;
27: Vec x,b;
28: DM da;
29: Mat A,Atrans;
30: PetscInt dof=1,M=-8;
31: PetscBool flg,trans=PETSC_FALSE;
33: PetscInitialize(&argc,&argv,(char *)0,help);
34: PetscOptionsGetInt(PETSC_NULL,"-dof",&dof,PETSC_NULL);
35: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
36: PetscOptionsGetBool(PETSC_NULL,"-trans",&trans,PETSC_NULL);
38: DMDACreate(PETSC_COMM_WORLD,&da);
39: DMDASetDim(da,3);
40: DMDASetBoundaryType(da,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE);
41: DMDASetStencilType(da,DMDA_STENCIL_STAR);
42: DMDASetSizes(da,M,M,M);
43: DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
44: DMDASetDof(da,dof);
45: DMDASetStencilWidth(da,1);
46: DMDASetOwnershipRanges(da,PETSC_NULL,PETSC_NULL,PETSC_NULL);
47: DMSetFromOptions(da);
48: DMSetUp(da);
50: DMCreateGlobalVector(da,&x);
51: DMCreateGlobalVector(da,&b);
52: ComputeRHS(da,b);
53: DMCreateMatrix(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: PetscOptionsGetBool(PETSC_NULL, "-test_sbaij1", &flg,PETSC_NULL);
66: if (flg){
67: Mat sA;
68: PetscBool issymm;
69: MatIsTranspose(A,A,0.0,&issymm);
70: if (issymm) {
71: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
72: } else {
73: printf("Warning: A is non-symmetric\n");
74: }
75: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
76: MatDestroy(&A);
77: A = sA;
78: }
80: KSPCreate(PETSC_COMM_WORLD,&ksp);
81: KSPSetFromOptions(ksp);
82: KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
83: KSPGetPC(ksp,&pc);
84: PCSetDM(pc,(DM)da);
85:
86: if (trans) {
87: KSPSolveTranspose(ksp,b,x);
88: } else {
89: KSPSolve(ksp,b,x);
90: }
92: /* check final residual */
93: flg = PETSC_FALSE;
94: PetscOptionsGetBool(PETSC_NULL, "-check_final_residual", &flg,PETSC_NULL);
95: if (flg){
96: Vec b1;
97: PetscReal norm;
98: KSPGetSolution(ksp,&x);
99: VecDuplicate(b,&b1);
100: MatMult(A,x,b1);
101: VecAXPY(b1,-1.0,b);
102: VecNorm(b1,NORM_2,&norm);
103: PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
104: VecDestroy(&b1);
105: }
106:
107: KSPDestroy(&ksp);
108: VecDestroy(&x);
109: VecDestroy(&b);
110: MatDestroy(&A);
111: DMDestroy(&da);
112: PetscFinalize();
113: return 0;
114: }
118: PetscErrorCode ComputeRHS(DM da,Vec b)
119: {
121: PetscInt mx,my,mz;
122: PetscScalar h;
125: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
126: h = 1.0/((mx-1)*(my-1)*(mz-1));
127: VecSet(b,h);
128: return(0);
129: }
130:
133: PetscErrorCode ComputeMatrix(DM da,Mat B)
134: {
136: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
137: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
138: MatStencil row,col;
139:
141: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
142: /* For simplicity, this example only works on mx=my=mz */
143: if ( mx != my || mx != mz) SETERRQ3(PETSC_COMM_SELF,1,"This example only works with mx %d = my %d = mz %d\n",mx,my,mz);
145: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
146: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
148: PetscMalloc((2*dof*dof+1)*sizeof(PetscScalar),&v);
149: v_neighbor = v + dof*dof;
150: PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
151: k3 = 0;
152: for (k1=0; k1<dof; k1++){
153: for (k2=0; k2<dof; k2++){
154: if (k1 == k2){
155: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
156: v_neighbor[k3] = -HxHydHz;
157: } else {
158: v[k3] = k1/(dof*dof); ;
159: v_neighbor[k3] = k2/(dof*dof);
160: }
161: k3++;
162: }
163: }
164: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
165:
166: for (k=zs; k<zs+zm; k++){
167: for (j=ys; j<ys+ym; j++){
168: for(i=xs; i<xs+xm; i++){
169: row.i = i; row.j = j; row.k = k;
170: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){ /* boudary points */
171: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
172: } else { /* interior points */
173: /* center */
174: col.i = i; col.j = j; col.k = k;
175: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
176:
177: /* x neighbors */
178: col.i = i-1; col.j = j; col.k = k;
179: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
180: col.i = i+1; col.j = j; col.k = k;
181: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
182:
183: /* y neighbors */
184: col.i = i; col.j = j-1; col.k = k;
185: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
186: col.i = i; col.j = j+1; col.k = k;
187: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
188:
189: /* z neighbors */
190: col.i = i; col.j = j; col.k = k-1;
191: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
192: col.i = i; col.j = j; col.k = k+1;
193: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
194: }
195: }
196: }
197: }
198: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
199: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
200: PetscFree(v);
201: return(0);
202: }