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: PetscFunctionBeginUser;
31: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
32: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
33: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
34: PetscCall(PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL));
36: PetscCall(DMDACreate(PETSC_COMM_WORLD, &da));
37: PetscCall(DMSetDimension(da, 3));
38: PetscCall(DMDASetBoundaryType(da, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE));
39: PetscCall(DMDASetStencilType(da, DMDA_STENCIL_STAR));
40: PetscCall(DMDASetSizes(da, M, M, M));
41: PetscCall(DMDASetNumProcs(da, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
42: PetscCall(DMDASetDof(da, dof));
43: PetscCall(DMDASetStencilWidth(da, 1));
44: PetscCall(DMDASetOwnershipRanges(da, NULL, NULL, NULL));
45: PetscCall(DMSetFromOptions(da));
46: PetscCall(DMSetUp(da));
48: PetscCall(DMCreateGlobalVector(da, &x));
49: PetscCall(DMCreateGlobalVector(da, &b));
50: PetscCall(ComputeRHS(da, b));
51: PetscCall(DMSetMatType(da, MATBAIJ));
52: PetscCall(DMSetFromOptions(da));
53: PetscCall(DMCreateMatrix(da, &A));
54: PetscCall(ComputeMatrix(da, A));
56: /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
57: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
58: PetscCall(MatAXPY(A, 1.0, Atrans, DIFFERENT_NONZERO_PATTERN));
59: PetscCall(MatScale(A, 0.5));
60: PetscCall(MatDestroy(&Atrans));
62: /* Test sbaij matrix */
63: flg = PETSC_FALSE;
64: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_sbaij1", &flg, NULL));
65: if (flg) {
66: Mat sA;
67: PetscBool issymm;
68: PetscCall(MatIsTranspose(A, A, 0.0, &issymm));
69: if (issymm) {
70: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
71: } else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: A is non-symmetric\n"));
72: PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));
73: PetscCall(MatDestroy(&A));
74: A = sA;
75: }
77: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
78: PetscCall(KSPSetFromOptions(ksp));
79: PetscCall(KSPSetOperators(ksp, A, A));
80: PetscCall(KSPGetPC(ksp, &pc));
81: PetscCall(PCSetDM(pc, (DM)da));
83: if (trans) {
84: PetscCall(KSPSolveTranspose(ksp, b, x));
85: } else {
86: PetscCall(KSPSolve(ksp, b, x));
87: }
89: /* check final residual */
90: flg = PETSC_FALSE;
91: PetscCall(PetscOptionsGetBool(NULL, NULL, "-check_final_residual", &flg, NULL));
92: if (flg) {
93: Vec b1;
94: PetscReal norm;
95: PetscCall(KSPGetSolution(ksp, &x));
96: PetscCall(VecDuplicate(b, &b1));
97: PetscCall(MatMult(A, x, b1));
98: PetscCall(VecAXPY(b1, -1.0, b));
99: PetscCall(VecNorm(b1, NORM_2, &norm));
100: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final residual %g\n", (double)norm));
101: PetscCall(VecDestroy(&b1));
102: }
104: PetscCall(KSPDestroy(&ksp));
105: PetscCall(VecDestroy(&x));
106: PetscCall(VecDestroy(&b));
107: PetscCall(MatDestroy(&A));
108: PetscCall(DMDestroy(&da));
109: PetscCall(PetscFinalize());
110: return 0;
111: }
113: PetscErrorCode ComputeRHS(DM da, Vec b)
114: {
115: PetscInt mx, my, mz;
116: PetscScalar h;
118: PetscFunctionBeginUser;
119: PetscCall(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: PetscCall(VecSet(b, h));
122: PetscFunctionReturn(PETSC_SUCCESS);
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;
131: PetscFunctionBeginUser;
132: PetscCall(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 */
134: PetscCheck(mx == my && mx == mz, PETSC_COMM_SELF, PETSC_ERR_SUP, "This example only works with mx %" PetscInt_FMT " = my %" PetscInt_FMT " = mz %" PetscInt_FMT, 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: PetscCall(PetscMalloc1(2 * dof * dof + 1, &v));
144: v_neighbor = v + dof * dof;
145: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
181: col.i = i + 1;
182: col.j = j;
183: col.k = k;
184: PetscCall(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: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
191: col.i = i;
192: col.j = j + 1;
193: col.k = k;
194: PetscCall(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: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
201: col.i = i;
202: col.j = j;
203: col.k = k + 1;
204: PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
205: }
206: }
207: }
208: }
209: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
210: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
211: PetscCall(PetscFree(v));
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*TEST
217: test:
218: args: -ksp_monitor_short -dm_mat_type sbaij -ksp_monitor_short -pc_type cholesky -ksp_view
220: TEST*/