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