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;
31: PetscInitialize(&argc, &argv, (char *)0, help);
32: PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL);
33: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
34: PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL);
36: DMDACreate(PETSC_COMM_WORLD, &da);
37: DMSetDimension(da, 3);
38: DMDASetBoundaryType(da, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE);
39: DMDASetStencilType(da, DMDA_STENCIL_STAR);
40: DMDASetSizes(da, M, M, M);
41: DMDASetNumProcs(da, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
42: DMDASetDof(da, dof);
43: DMDASetStencilWidth(da, 1);
44: DMDASetOwnershipRanges(da, NULL, NULL, NULL);
45: DMSetFromOptions(da);
46: DMSetUp(da);
48: DMCreateGlobalVector(da, &x);
49: DMCreateGlobalVector(da, &b);
50: ComputeRHS(da, b);
51: DMSetMatType(da, MATBAIJ);
52: DMSetFromOptions(da);
53: DMCreateMatrix(da, &A);
54: ComputeMatrix(da, A);
56: /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
57: MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans);
58: MatAXPY(A, 1.0, Atrans, DIFFERENT_NONZERO_PATTERN);
59: MatScale(A, 0.5);
60: MatDestroy(&Atrans);
62: /* Test sbaij matrix */
63: flg = PETSC_FALSE;
64: PetscOptionsGetBool(NULL, NULL, "-test_sbaij1", &flg, NULL);
65: if (flg) {
66: Mat sA;
67: PetscBool issymm;
68: MatIsTranspose(A, A, 0.0, &issymm);
69: if (issymm) {
70: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
71: } else PetscPrintf(PETSC_COMM_WORLD, "Warning: A is non-symmetric\n");
72: MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA);
73: MatDestroy(&A);
74: A = sA;
75: }
77: KSPCreate(PETSC_COMM_WORLD, &ksp);
78: KSPSetFromOptions(ksp);
79: KSPSetOperators(ksp, A, A);
80: KSPGetPC(ksp, &pc);
81: PCSetDM(pc, (DM)da);
83: if (trans) {
84: KSPSolveTranspose(ksp, b, x);
85: } else {
86: KSPSolve(ksp, b, x);
87: }
89: /* check final residual */
90: flg = PETSC_FALSE;
91: PetscOptionsGetBool(NULL, NULL, "-check_final_residual", &flg, NULL);
92: if (flg) {
93: Vec b1;
94: PetscReal norm;
95: KSPGetSolution(ksp, &x);
96: VecDuplicate(b, &b1);
97: MatMult(A, x, b1);
98: VecAXPY(b1, -1.0, b);
99: VecNorm(b1, NORM_2, &norm);
100: PetscPrintf(PETSC_COMM_WORLD, "Final residual %g\n", (double)norm);
101: VecDestroy(&b1);
102: }
104: KSPDestroy(&ksp);
105: VecDestroy(&x);
106: VecDestroy(&b);
107: MatDestroy(&A);
108: DMDestroy(&da);
109: PetscFinalize();
110: return 0;
111: }
113: PetscErrorCode ComputeRHS(DM da, Vec b)
114: {
115: PetscInt mx, my, mz;
116: PetscScalar h;
119: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
120: h = 1.0 / ((mx - 1) * (my - 1) * (mz - 1));
121: VecSet(b, h);
122: return 0;
123: }
125: PetscErrorCode ComputeMatrix(DM da, Mat B)
126: {
127: PetscInt i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs, dof, k1, k2, k3;
128: PetscScalar *v, *v_neighbor, Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy;
129: MatStencil row, col;
132: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
133: /* For simplicity, this example only works on mx=my=mz */
136: Hx = 1.0 / (PetscReal)(mx - 1);
137: Hy = 1.0 / (PetscReal)(my - 1);
138: Hz = 1.0 / (PetscReal)(mz - 1);
139: HxHydHz = Hx * Hy / Hz;
140: HxHzdHy = Hx * Hz / Hy;
141: HyHzdHx = Hy * Hz / Hx;
143: PetscMalloc1(2 * dof * dof + 1, &v);
144: v_neighbor = v + dof * dof;
145: PetscArrayzero(v, 2 * dof * dof + 1);
146: k3 = 0;
147: for (k1 = 0; k1 < dof; k1++) {
148: for (k2 = 0; k2 < dof; k2++) {
149: if (k1 == k2) {
150: v[k3] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx);
151: v_neighbor[k3] = -HxHydHz;
152: } else {
153: v[k3] = k1 / (dof * dof);
154: v_neighbor[k3] = k2 / (dof * dof);
155: }
156: k3++;
157: }
158: }
159: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
161: for (k = zs; k < zs + zm; k++) {
162: for (j = ys; j < ys + ym; j++) {
163: for (i = xs; i < xs + xm; i++) {
164: row.i = i;
165: row.j = j;
166: row.k = k;
167: if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || k == mz - 1) { /* boundary points */
168: MatSetValuesBlockedStencil(B, 1, &row, 1, &row, v, INSERT_VALUES);
169: } else { /* interior points */
170: /* center */
171: col.i = i;
172: col.j = j;
173: col.k = k;
174: MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v, INSERT_VALUES);
176: /* x neighbors */
177: col.i = i - 1;
178: col.j = j;
179: col.k = k;
180: MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES);
181: col.i = i + 1;
182: col.j = j;
183: col.k = k;
184: MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES);
186: /* y neighbors */
187: col.i = i;
188: col.j = j - 1;
189: col.k = k;
190: MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES);
191: col.i = i;
192: col.j = j + 1;
193: col.k = k;
194: MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES);
196: /* z neighbors */
197: col.i = i;
198: col.j = j;
199: col.k = k - 1;
200: MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES);
201: col.i = i;
202: col.j = j;
203: col.k = k + 1;
204: MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES);
205: }
206: }
207: }
208: }
209: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
210: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
211: PetscFree(v);
212: return 0;
213: }
215: /*TEST
217: test:
218: args: -ksp_monitor_short -dm_mat_type sbaij -ksp_monitor_short -pc_type cholesky -ksp_view
220: TEST*/