Actual source code: ex32.c
petsc-3.5.4 2015-05-23
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);
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(NULL,"-dof",&dof,NULL);
35: PetscOptionsGetInt(NULL,"-M",&M,NULL);
36: PetscOptionsGetBool(NULL,"-trans",&trans,NULL);
38: DMDACreate(PETSC_COMM_WORLD,&da);
39: DMDASetDim(da,3);
40: DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_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,NULL,NULL,NULL);
47: DMSetFromOptions(da);
48: DMSetUp(da);
50: DMCreateGlobalVector(da,&x);
51: DMCreateGlobalVector(da,&b);
52: ComputeRHS(da,b);
53: DMSetMatType(da,MATBAIJ);
54: DMCreateMatrix(da,&A);
55: ComputeMatrix(da,A);
58: /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
59: MatTranspose(A,MAT_INITIAL_MATRIX,&Atrans);
60: MatAXPY(A,1.0,Atrans,DIFFERENT_NONZERO_PATTERN);
61: MatScale(A,0.5);
62: MatDestroy(&Atrans);
64: /* Test sbaij matrix */
65: flg = PETSC_FALSE;
66: PetscOptionsGetBool(NULL, "-test_sbaij1", &flg,NULL);
67: if (flg) {
68: Mat sA;
69: PetscBool issymm;
70: MatIsTranspose(A,A,0.0,&issymm);
71: if (issymm) {
72: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
73: } else printf("Warning: A is non-symmetric\n");
74: MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
75: MatDestroy(&A);
76: A = sA;
77: }
79: KSPCreate(PETSC_COMM_WORLD,&ksp);
80: KSPSetFromOptions(ksp);
81: KSPSetOperators(ksp,A,A);
82: KSPGetPC(ksp,&pc);
83: PCSetDM(pc,(DM)da);
85: if (trans) {
86: KSPSolveTranspose(ksp,b,x);
87: } else {
88: KSPSolve(ksp,b,x);
89: }
91: /* check final residual */
92: flg = PETSC_FALSE;
93: PetscOptionsGetBool(NULL, "-check_final_residual", &flg,NULL);
94: if (flg) {
95: Vec b1;
96: PetscReal norm;
97: KSPGetSolution(ksp,&x);
98: VecDuplicate(b,&b1);
99: MatMult(A,x,b1);
100: VecAXPY(b1,-1.0,b);
101: VecNorm(b1,NORM_2,&norm);
102: PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
103: VecDestroy(&b1);
104: }
106: KSPDestroy(&ksp);
107: VecDestroy(&x);
108: VecDestroy(&b);
109: MatDestroy(&A);
110: DMDestroy(&da);
111: PetscFinalize();
112: return 0;
113: }
117: PetscErrorCode ComputeRHS(DM da,Vec b)
118: {
120: PetscInt mx,my,mz;
121: PetscScalar h;
124: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
125: h = 1.0/((mx-1)*(my-1)*(mz-1));
126: VecSet(b,h);
127: return(0);
128: }
132: PetscErrorCode ComputeMatrix(DM da,Mat B)
133: {
135: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,dof,k1,k2,k3;
136: PetscScalar *v,*v_neighbor,Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
137: MatStencil row,col;
140: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
141: /* For simplicity, this example only works on mx=my=mz */
142: 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);
144: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
145: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
147: PetscMalloc1((2*dof*dof+1),&v);
148: v_neighbor = v + dof*dof;
149: PetscMemzero(v,(2*dof*dof+1)*sizeof(PetscScalar));
150: k3 = 0;
151: for (k1=0; k1<dof; k1++) {
152: for (k2=0; k2<dof; k2++) {
153: if (k1 == k2) {
154: v[k3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
155: v_neighbor[k3] = -HxHydHz;
156: } else {
157: v[k3] = k1/(dof*dof);;
158: v_neighbor[k3] = k2/(dof*dof);
159: }
160: k3++;
161: }
162: }
163: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
165: for (k=zs; k<zs+zm; k++) {
166: for (j=ys; j<ys+ym; j++) {
167: for (i=xs; i<xs+xm; i++) {
168: row.i = i; row.j = j; row.k = k;
169: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) { /* boudary points */
170: MatSetValuesBlockedStencil(B,1,&row,1,&row,v,INSERT_VALUES);
171: } else { /* interior points */
172: /* center */
173: col.i = i; col.j = j; col.k = k;
174: MatSetValuesBlockedStencil(B,1,&row,1,&col,v,INSERT_VALUES);
176: /* x neighbors */
177: col.i = i-1; col.j = j; col.k = k;
178: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
179: col.i = i+1; col.j = j; col.k = k;
180: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
182: /* y neighbors */
183: col.i = i; col.j = j-1; col.k = k;
184: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
185: col.i = i; col.j = j+1; col.k = k;
186: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
188: /* z neighbors */
189: col.i = i; col.j = j; col.k = k-1;
190: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
191: col.i = i; col.j = j; col.k = k+1;
192: MatSetValuesBlockedStencil(B,1,&row,1,&col,v_neighbor,INSERT_VALUES);
193: }
194: }
195: }
196: }
197: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
198: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
199: PetscFree(v);
200: return(0);
201: }