Actual source code: ex32.c
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 <petscdm.h>
14: #include <petscdmda.h>
15: #include <petscksp.h>
17: extern PetscErrorCode ComputeMatrix(DM,Mat);
18: extern PetscErrorCode ComputeRHS(DM,Vec);
20: int main(int argc,char **argv)
21: {
22: KSP ksp;
23: PC pc;
24: Vec x,b;
25: DM da;
26: Mat A,Atrans;
27: PetscInt dof=1,M=8;
28: PetscBool flg,trans=PETSC_FALSE;
30: PetscInitialize(&argc,&argv,(char*)0,help);
31: PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL);
32: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
33: PetscOptionsGetBool(NULL,NULL,"-trans",&trans,NULL);
35: DMDACreate(PETSC_COMM_WORLD,&da);
36: DMSetDimension(da,3);
37: DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE);
38: DMDASetStencilType(da,DMDA_STENCIL_STAR);
39: DMDASetSizes(da,M,M,M);
40: DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
41: DMDASetDof(da,dof);
42: DMDASetStencilWidth(da,1);
43: DMDASetOwnershipRanges(da,NULL,NULL,NULL);
44: DMSetFromOptions(da);
45: DMSetUp(da);
47: DMCreateGlobalVector(da,&x);
48: DMCreateGlobalVector(da,&b);
49: ComputeRHS(da,b);
50: DMSetMatType(da,MATBAIJ);
51: DMSetFromOptions(da);
52: DMCreateMatrix(da,&A);
53: ComputeMatrix(da,A);
55: /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
56: MatTranspose(A,MAT_INITIAL_MATRIX,&Atrans);
57: MatAXPY(A,1.0,Atrans,DIFFERENT_NONZERO_PATTERN);
58: MatScale(A,0.5);
59: MatDestroy(&Atrans);
61: /* Test sbaij matrix */
62: flg = PETSC_FALSE;
63: PetscOptionsGetBool(NULL,NULL, "-test_sbaij1", &flg,NULL);
64: if (flg) {
65: Mat sA;
66: PetscBool issymm;
67: MatIsTranspose(A,A,0.0,&issymm);
68: if (issymm) {
69: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
70: } else PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric\n");
71: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
72: MatDestroy(&A);
73: A = sA;
74: }
76: KSPCreate(PETSC_COMM_WORLD,&ksp);
77: KSPSetFromOptions(ksp);
78: KSPSetOperators(ksp,A,A);
79: KSPGetPC(ksp,&pc);
80: PCSetDM(pc,(DM)da);
82: if (trans) {
83: KSPSolveTranspose(ksp,b,x);
84: } else {
85: KSPSolve(ksp,b,x);
86: }
88: /* check final residual */
89: flg = PETSC_FALSE;
90: PetscOptionsGetBool(NULL,NULL, "-check_final_residual", &flg,NULL);
91: if (flg) {
92: Vec b1;
93: PetscReal norm;
94: KSPGetSolution(ksp,&x);
95: VecDuplicate(b,&b1);
96: MatMult(A,x,b1);
97: VecAXPY(b1,-1.0,b);
98: VecNorm(b1,NORM_2,&norm);
99: PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
100: VecDestroy(&b1);
101: }
103: KSPDestroy(&ksp);
104: VecDestroy(&x);
105: VecDestroy(&b);
106: MatDestroy(&A);
107: DMDestroy(&da);
108: PetscFinalize();
109: return 0;
110: }
112: PetscErrorCode ComputeRHS(DM da,Vec b)
113: {
114: PetscInt mx,my,mz;
115: PetscScalar h;
118: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,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: }
124: PetscErrorCode ComputeMatrix(DM da,Mat B)
125: {
126: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
127: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
128: MatStencil row,col;
131: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
132: /* For simplicity, this example only works on mx=my=mz */
135: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
136: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
138: PetscMalloc1(2*dof*dof+1,&v);
139: v_neighbor = v + dof*dof;
140: PetscArrayzero(v,2*dof*dof+1);
141: k3 = 0;
142: for (k1=0; k1<dof; k1++) {
143: for (k2=0; k2<dof; k2++) {
144: if (k1 == k2) {
145: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
146: v_neighbor[k3] = -HxHydHz;
147: } else {
148: v[k3] = k1/(dof*dof);
149: v_neighbor[k3] = k2/(dof*dof);
150: }
151: k3++;
152: }
153: }
154: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
156: for (k=zs; k<zs+zm; k++) {
157: for (j=ys; j<ys+ym; j++) {
158: for (i=xs; i<xs+xm; i++) {
159: row.i = i; row.j = j; row.k = k;
160: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { /* boundary points */
161: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
162: } else { /* interior points */
163: /* center */
164: col.i = i; col.j = j; col.k = k;
165: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
167: /* x neighbors */
168: col.i = i-1; col.j = j; col.k = k;
169: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
170: col.i = i+1; col.j = j; col.k = k;
171: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
173: /* y neighbors */
174: col.i = i; col.j = j-1; col.k = k;
175: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
176: col.i = i; col.j = j+1; col.k = k;
177: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
179: /* z neighbors */
180: col.i = i; col.j = j; col.k = k-1;
181: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
182: col.i = i; col.j = j; col.k = k+1;
183: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
184: }
185: }
186: }
187: }
188: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
189: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
190: PetscFree(v);
191: return 0;
192: }
194: /*TEST
196: test:
197: args: -ksp_monitor_short -dm_mat_type sbaij -ksp_monitor_short -pc_type cholesky -ksp_view
199: TEST*/