Actual source code: ex26.c
1: static char help[] = "Solves Laplacian with multigrid. Tests block API for PCMG\n\
2: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
3: -my <yg>, where <yg> = number of grid points in the y-direction\n\
4: -Nx <npx>, where <npx> = number of processors in the x-direction\n\
5: -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";
7: /* Modified from ~src/ksp/tests/ex19.c. Used for testing ML 6.2 interface.
9: This problem is modeled by
10: the partial differential equation
12: -Laplacian u = g, 0 < x,y < 1,
14: with boundary conditions
16: u = 0 for x = 0, x = 1, y = 0, y = 1.
18: A finite difference approximation with the usual 5-point stencil
19: is used to discretize the boundary value problem to obtain a linear
20: system of equations.
22: Usage: ./ex26 -ksp_monitor -pc_type ml
23: -mg_coarse_ksp_max_it 10
24: -mg_levels_1_ksp_max_it 10 -mg_levels_2_ksp_max_it 10
25: -mg_fine_ksp_max_it 10
26: */
28: #include <petscksp.h>
29: #include <petscdm.h>
30: #include <petscdmda.h>
32: /* User-defined application contexts */
33: typedef struct {
34: PetscInt mx, my; /* number grid points in x and y direction */
35: Vec localX, localF; /* local vectors with ghost region */
36: DM da;
37: Vec x, b, r; /* global vectors */
38: Mat J; /* Jacobian on grid */
39: Mat A, P, R;
40: KSP ksp;
41: } GridCtx;
43: static PetscErrorCode FormJacobian_Grid(GridCtx *, Mat);
45: int main(int argc, char **argv)
46: {
47: PetscInt i, its, Nx = PETSC_DECIDE, Ny = PETSC_DECIDE, nlocal, nrhs = 1;
48: PetscScalar one = 1.0;
49: Mat A, B, X;
50: GridCtx fine_ctx;
51: KSP ksp;
52: PetscBool Brand = PETSC_FALSE, transpose = PETSC_FALSE, flg;
54: PetscFunctionBeginUser;
55: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
56: /* set up discretization matrix for fine grid */
57: fine_ctx.mx = 9;
58: fine_ctx.my = 9;
59: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &fine_ctx.mx, NULL));
60: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &fine_ctx.my, NULL));
61: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
62: PetscCall(PetscOptionsGetInt(NULL, NULL, "-Nx", &Nx, NULL));
63: PetscCall(PetscOptionsGetInt(NULL, NULL, "-Ny", &Ny, NULL));
64: PetscCall(PetscOptionsGetBool(NULL, NULL, "-rand", &Brand, NULL));
65: PetscCall(PetscOptionsGetBool(NULL, NULL, "-transpose", &transpose, NULL));
66: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Fine grid size %" PetscInt_FMT " by %" PetscInt_FMT "\n", fine_ctx.mx, fine_ctx.my));
68: /* Set up distributed array for fine grid */
69: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, fine_ctx.mx, fine_ctx.my, Nx, Ny, 1, 1, NULL, NULL, &fine_ctx.da));
70: PetscCall(DMSetFromOptions(fine_ctx.da));
71: PetscCall(DMSetUp(fine_ctx.da));
72: PetscCall(DMCreateGlobalVector(fine_ctx.da, &fine_ctx.x));
73: PetscCall(VecDuplicate(fine_ctx.x, &fine_ctx.b));
74: PetscCall(VecGetLocalSize(fine_ctx.x, &nlocal));
75: PetscCall(DMCreateLocalVector(fine_ctx.da, &fine_ctx.localX));
76: PetscCall(VecDuplicate(fine_ctx.localX, &fine_ctx.localF));
77: PetscCall(DMCreateMatrix(fine_ctx.da, &A));
78: PetscCall(FormJacobian_Grid(&fine_ctx, A));
80: /* create linear solver */
81: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
82: PetscCall(KSPSetDM(ksp, fine_ctx.da));
83: PetscCall(KSPSetDMActive(ksp, KSP_DMACTIVE_ALL, PETSC_FALSE));
85: /* set values for rhs vector */
86: PetscCall(VecSet(fine_ctx.b, one));
88: /* set options, then solve system */
89: PetscCall(KSPSetFromOptions(ksp)); /* calls PCSetFromOptions_ML if 'pc_type=ml' */
90: PetscCall(KSPSetOperators(ksp, A, A));
91: PetscCall(KSPSolve(ksp, fine_ctx.b, fine_ctx.x));
92: PetscCall(VecViewFromOptions(fine_ctx.x, NULL, "-debug"));
93: PetscCall(KSPGetIterationNumber(ksp, &its));
94: PetscCall(KSPGetIterationNumber(ksp, &its));
95: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %" PetscInt_FMT "\n", its));
97: /* test multiple right-hand side */
98: PetscCall(MatCreateDense(PETSC_COMM_WORLD, nlocal, PETSC_DECIDE, fine_ctx.mx * fine_ctx.my, nrhs, NULL, &B));
99: PetscCall(MatSetOptionsPrefix(B, "rhs_"));
100: PetscCall(MatSetFromOptions(B));
101: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &X));
102: /* start the -ksp_initial_guess_nonzero runs from a known state, MatDuplicate() zeroing is not guaranteed for all -rhs_mat_type values */
103: PetscCall(MatZeroEntries(X));
104: if (Brand) {
105: PetscCall(MatSetRandom(B, NULL));
106: } else {
107: PetscScalar *b;
109: PetscCall(MatDenseGetArrayWrite(B, &b));
110: for (i = 0; i < nlocal * nrhs; i++) b[i] = 1.0;
111: PetscCall(MatDenseRestoreArrayWrite(B, &b));
112: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
113: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
114: }
115: if (transpose) PetscCall(KSPMatSolveTranspose(ksp, B, X));
116: else PetscCall(KSPMatSolve(ksp, B, X));
117: PetscCall(MatViewFromOptions(X, NULL, "-debug"));
119: PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPPREONLY, &flg));
120: if ((flg || nrhs == 1) && !Brand && !transpose) { /* the operator is not symmetric, so KSPMatSolveTranspose() does not return the solution of KSPSolve() */
121: const PetscScalar *xx, *XX;
123: PetscCall(VecGetArrayRead(fine_ctx.x, &xx));
124: PetscCall(MatDenseGetArrayRead(X, &XX));
125: for (PetscInt n = 0; n < nrhs; n++) {
126: for (i = 0; i < nlocal; i++) {
127: if (PetscAbsScalar(xx[i] - XX[nlocal * n + i]) > PETSC_SMALL) {
128: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Error local solve %" PetscInt_FMT ", entry %" PetscInt_FMT " -> %g + i %g != %g + i %g\n", PetscGlobalRank, n, i, (double)PetscRealPart(xx[i]), (double)PetscImaginaryPart(xx[i]), (double)PetscRealPart(XX[i]), (double)PetscImaginaryPart(XX[i])));
129: }
130: }
131: }
132: PetscCall(MatDenseRestoreArrayRead(X, &XX));
133: PetscCall(VecRestoreArrayRead(fine_ctx.x, &xx));
134: }
136: /* free data structures */
137: PetscCall(VecDestroy(&fine_ctx.x));
138: PetscCall(VecDestroy(&fine_ctx.b));
139: PetscCall(DMDestroy(&fine_ctx.da));
140: PetscCall(VecDestroy(&fine_ctx.localX));
141: PetscCall(VecDestroy(&fine_ctx.localF));
142: PetscCall(MatDestroy(&A));
143: PetscCall(MatDestroy(&B));
144: PetscCall(MatDestroy(&X));
145: PetscCall(KSPDestroy(&ksp));
147: PetscCall(PetscFinalize());
148: return 0;
149: }
151: PetscErrorCode FormJacobian_Grid(GridCtx *grid, Mat jac)
152: {
153: PetscInt i, j, row, mx, my, xs, ys, xm, ym, Xs, Ys, Xm, Ym, col[5];
154: PetscInt grow;
155: const PetscInt *ltog;
156: PetscScalar two = 2.0, one = 1.0, v[5], hx, hy, hxdhy, hydhx, value;
157: ISLocalToGlobalMapping ltogm;
159: PetscFunctionBeginUser;
160: mx = grid->mx;
161: my = grid->my;
162: hx = one / (PetscReal)(mx - 1);
163: hy = one / (PetscReal)(my - 1);
164: hxdhy = hx / hy;
165: hydhx = hy / hx;
167: /* Get ghost points */
168: PetscCall(DMDAGetCorners(grid->da, &xs, &ys, 0, &xm, &ym, 0));
169: PetscCall(DMDAGetGhostCorners(grid->da, &Xs, &Ys, 0, &Xm, &Ym, 0));
170: PetscCall(DMGetLocalToGlobalMapping(grid->da, <ogm));
171: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, <og));
173: /* Evaluate Jacobian of function */
174: for (j = ys; j < ys + ym; j++) {
175: row = (j - Ys) * Xm + xs - Xs - 1;
176: for (i = xs; i < xs + xm; i++) {
177: row++;
178: grow = ltog[row];
179: if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
180: v[0] = -hxdhy;
181: col[0] = ltog[row - Xm];
182: v[1] = -hydhx;
183: col[1] = ltog[row - 1];
184: v[2] = two * (hydhx + hxdhy);
185: col[2] = grow;
186: v[3] = -hydhx;
187: col[3] = ltog[row + 1];
188: v[4] = -hxdhy;
189: col[4] = ltog[row + Xm];
190: PetscCall(MatSetValues(jac, 1, &grow, 5, col, v, INSERT_VALUES));
191: } else if ((i > 0 && i < mx - 1) || (j > 0 && j < my - 1)) {
192: value = .5 * two * (hydhx + hxdhy);
193: PetscCall(MatSetValues(jac, 1, &grow, 1, &grow, &value, INSERT_VALUES));
194: } else {
195: value = .25 * two * (hydhx + hxdhy);
196: PetscCall(MatSetValues(jac, 1, &grow, 1, &grow, &value, INSERT_VALUES));
197: }
198: }
199: }
200: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, <og));
201: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
202: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*TEST
208: test:
209: args: -ksp_monitor
211: test:
212: suffix: 2
213: args: -ksp_monitor
214: nsize: 3
216: test:
217: suffix: ml_1
218: args: -ksp_monitor -pc_type ml -mat_no_inode
219: nsize: 3
220: requires: ml
222: test:
223: suffix: ml_2
224: args: -ksp_monitor -pc_type ml -mat_no_inode -ksp_max_it 3
225: nsize: 3
226: requires: ml
228: test:
229: suffix: ml_3
230: args: -ksp_monitor -pc_type ml -mat_no_inode -pc_mg_type ADDITIVE -ksp_max_it 7
231: nsize: 1
232: requires: ml
234: test:
235: suffix: cycles
236: nsize: {{1 2}}
237: args: -ksp_view_final_residual -pc_type mg -mx 5 -my 5 -pc_mg_levels 3 -pc_mg_galerkin -ksp_monitor -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -pc_mg_type {{additive multiplicative full kaskade}separate output} -nrhs 1
239: test:
240: suffix: matcycles
241: nsize: {{1 2}}
242: args: -ksp_view_final_residual -ksp_type preonly -pc_type mg -mx 5 -my 5 -pc_mg_levels 3 -pc_mg_galerkin -ksp_monitor -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -pc_mg_type {{additive multiplicative full kaskade}separate output} -nrhs 7 -ksp_matsolve_batch_size {{4 7}separate output}
244: test:
245: suffix: matcycles_richardson
246: nsize: {{1 2}}
247: args: -ksp_view_final_residual -ksp_type richardson -pc_type mg -mx 5 -my 5 -pc_mg_levels 3 -pc_mg_galerkin -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -nrhs 7 -ksp_matsolve_batch_size {{4 7}separate output}
249: test:
250: suffix: matcycles_richardson_guess
251: nsize: {{1 2}}
252: args: -ksp_view_final_residual -ksp_type richardson -ksp_initial_guess_nonzero -ksp_norm_type unpreconditioned -pc_type mg -mx 5 -my 5 -pc_mg_levels 3 -pc_mg_galerkin -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -nrhs 7 -ksp_matsolve_batch_size {{4 7}separate output}
254: test:
255: suffix: matcycles_richardson_transpose
256: nsize: {{1 2}}
257: args: -ksp_view_final_residual -ksp_type richardson -transpose -pc_type jacobi -mx 5 -my 5 -nrhs 7 -ksp_matsolve_batch_size {{4 7}separate output}
259: test:
260: requires: ml
261: suffix: matcycles_ml
262: nsize: {{1 2}}
263: args: -ksp_view_final_residual -ksp_type preonly -pc_type ml -mx 5 -my 5 -ksp_monitor -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -pc_mg_type {{additive multiplicative full kaskade}separate output} -nrhs 7 -ksp_matsolve_batch_size {{4 7}separate output}
265: testset:
266: requires: hpddm
267: args: -ksp_view_final_residual -ksp_type hpddm -pc_type mg -pc_mg_levels 3 -pc_mg_galerkin -mx 5 -my 5 -ksp_monitor -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -nrhs 7
268: test:
269: suffix: matcycles_hpddm_mg
270: nsize: {{1 2}}
271: args: -pc_mg_type {{additive multiplicative full kaskade}separate output} -ksp_matsolve_batch_size {{4 7}separate output}
272: test:
273: requires: !__float128 !__fp16
274: suffix: hpddm_mg_mixed_precision
275: nsize: 2
276: output_file: output/ex26_matcycles_hpddm_mg_pc_mg_type-multiplicative_ksp_matsolve_batch_size-4.out
277: args: -ksp_matsolve_batch_size 4 -ksp_hpddm_precision {{single double}shared output}
278: test:
279: requires: __float128
280: suffix: hpddm_mg_mixed_precision___float128
281: nsize: 2
282: output_file: output/ex26_matcycles_hpddm_mg_pc_mg_type-multiplicative_ksp_matsolve_batch_size-4.out
283: args: -ksp_matsolve_batch_size 4 -ksp_hpddm_precision {{double __float128}shared output}
284: test:
285: requires: double defined(PETSC_HAVE_F2CBLASLAPACK___FLOAT128_BINDINGS)
286: suffix: hpddm_mg_mixed_precision_double
287: nsize: 2
288: output_file: output/ex26_matcycles_hpddm_mg_pc_mg_type-multiplicative_ksp_matsolve_batch_size-4.out
289: args: -ksp_matsolve_batch_size 4 -ksp_hpddm_precision __float128
290: test:
291: requires: single defined(PETSC_HAVE_F2CBLASLAPACK___FP16_BINDINGS)
292: suffix: hpddm_mg_mixed_precision_single
293: nsize: 2
294: output_file: output/ex26_matcycles_hpddm_mg_pc_mg_type-multiplicative_ksp_matsolve_batch_size-4.out
295: args: -ksp_matsolve_batch_size 4 -ksp_hpddm_precision __fp16 -ksp_rtol 1e-3
297: test:
298: requires: hpddm
299: nsize: {{1 2}}
300: suffix: matcycles_hpddm_ilu
301: args: -ksp_view_final_residual -ksp_type hpddm -pc_type redundant -redundant_pc_type ilu -mx 5 -my 5 -ksp_monitor -nrhs 7 -ksp_matsolve_batch_size {{4 7}separate output}
303: TEST*/