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_short -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, 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(PetscPrintf(PETSC_COMM_WORLD, "Fine grid size %" PetscInt_FMT " by %" PetscInt_FMT "\n", fine_ctx.mx, fine_ctx.my));

 67:   /* Set up distributed array for fine grid */
 68:   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));
 69:   PetscCall(DMSetFromOptions(fine_ctx.da));
 70:   PetscCall(DMSetUp(fine_ctx.da));
 71:   PetscCall(DMCreateGlobalVector(fine_ctx.da, &fine_ctx.x));
 72:   PetscCall(VecDuplicate(fine_ctx.x, &fine_ctx.b));
 73:   PetscCall(VecGetLocalSize(fine_ctx.x, &nlocal));
 74:   PetscCall(DMCreateLocalVector(fine_ctx.da, &fine_ctx.localX));
 75:   PetscCall(VecDuplicate(fine_ctx.localX, &fine_ctx.localF));
 76:   PetscCall(DMCreateMatrix(fine_ctx.da, &A));
 77:   PetscCall(FormJacobian_Grid(&fine_ctx, A));

 79:   /* create linear solver */
 80:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 81:   PetscCall(KSPSetDM(ksp, fine_ctx.da));
 82:   PetscCall(KSPSetDMActive(ksp, PETSC_FALSE));

 84:   /* set values for rhs vector */
 85:   PetscCall(VecSet(fine_ctx.b, one));

 87:   /* set options, then solve system */
 88:   PetscCall(KSPSetFromOptions(ksp)); /* calls PCSetFromOptions_ML if 'pc_type=ml' */
 89:   PetscCall(KSPSetOperators(ksp, A, A));
 90:   PetscCall(KSPSolve(ksp, fine_ctx.b, fine_ctx.x));
 91:   PetscCall(VecViewFromOptions(fine_ctx.x, NULL, "-debug"));
 92:   PetscCall(KSPGetIterationNumber(ksp, &its));
 93:   PetscCall(KSPGetIterationNumber(ksp, &its));
 94:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %" PetscInt_FMT "\n", its));

 96:   /* test multiple right-hand side */
 97:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, nlocal, PETSC_DECIDE, fine_ctx.mx * fine_ctx.my, nrhs, NULL, &B));
 98:   PetscCall(MatSetOptionsPrefix(B, "rhs_"));
 99:   PetscCall(MatSetFromOptions(B));
100:   PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &X));
101:   if (Brand) {
102:     PetscCall(MatSetRandom(B, NULL));
103:   } else {
104:     PetscScalar *b;

106:     PetscCall(MatDenseGetArrayWrite(B, &b));
107:     for (i = 0; i < nlocal * nrhs; i++) b[i] = 1.0;
108:     PetscCall(MatDenseRestoreArrayWrite(B, &b));
109:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
110:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
111:   }
112:   PetscCall(KSPMatSolve(ksp, B, X));
113:   PetscCall(MatViewFromOptions(X, NULL, "-debug"));

115:   PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPPREONLY, &flg));
116:   if ((flg || nrhs == 1) && !Brand) {
117:     PetscInt           n;
118:     const PetscScalar *xx, *XX;

120:     PetscCall(VecGetArrayRead(fine_ctx.x, &xx));
121:     PetscCall(MatDenseGetArrayRead(X, &XX));
122:     for (n = 0; n < nrhs; n++) {
123:       for (i = 0; i < nlocal; i++) {
124:         if (PetscAbsScalar(xx[i] - XX[nlocal * n + i]) > PETSC_SMALL) {
125:           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])));
126:         }
127:       }
128:     }
129:     PetscCall(MatDenseRestoreArrayRead(X, &XX));
130:     PetscCall(VecRestoreArrayRead(fine_ctx.x, &xx));
131:   }

133:   /* free data structures */
134:   PetscCall(VecDestroy(&fine_ctx.x));
135:   PetscCall(VecDestroy(&fine_ctx.b));
136:   PetscCall(DMDestroy(&fine_ctx.da));
137:   PetscCall(VecDestroy(&fine_ctx.localX));
138:   PetscCall(VecDestroy(&fine_ctx.localF));
139:   PetscCall(MatDestroy(&A));
140:   PetscCall(MatDestroy(&B));
141:   PetscCall(MatDestroy(&X));
142:   PetscCall(KSPDestroy(&ksp));

144:   PetscCall(PetscFinalize());
145:   return 0;
146: }

148: PetscErrorCode FormJacobian_Grid(GridCtx *grid, Mat jac)
149: {
150:   PetscInt               i, j, row, mx, my, xs, ys, xm, ym, Xs, Ys, Xm, Ym, col[5];
151:   PetscInt               grow;
152:   const PetscInt        *ltog;
153:   PetscScalar            two = 2.0, one = 1.0, v[5], hx, hy, hxdhy, hydhx, value;
154:   ISLocalToGlobalMapping ltogm;

156:   PetscFunctionBeginUser;
157:   mx    = grid->mx;
158:   my    = grid->my;
159:   hx    = one / (PetscReal)(mx - 1);
160:   hy    = one / (PetscReal)(my - 1);
161:   hxdhy = hx / hy;
162:   hydhx = hy / hx;

164:   /* Get ghost points */
165:   PetscCall(DMDAGetCorners(grid->da, &xs, &ys, 0, &xm, &ym, 0));
166:   PetscCall(DMDAGetGhostCorners(grid->da, &Xs, &Ys, 0, &Xm, &Ym, 0));
167:   PetscCall(DMGetLocalToGlobalMapping(grid->da, &ltogm));
168:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &ltog));

170:   /* Evaluate Jacobian of function */
171:   for (j = ys; j < ys + ym; j++) {
172:     row = (j - Ys) * Xm + xs - Xs - 1;
173:     for (i = xs; i < xs + xm; i++) {
174:       row++;
175:       grow = ltog[row];
176:       if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
177:         v[0]   = -hxdhy;
178:         col[0] = ltog[row - Xm];
179:         v[1]   = -hydhx;
180:         col[1] = ltog[row - 1];
181:         v[2]   = two * (hydhx + hxdhy);
182:         col[2] = grow;
183:         v[3]   = -hydhx;
184:         col[3] = ltog[row + 1];
185:         v[4]   = -hxdhy;
186:         col[4] = ltog[row + Xm];
187:         PetscCall(MatSetValues(jac, 1, &grow, 5, col, v, INSERT_VALUES));
188:       } else if ((i > 0 && i < mx - 1) || (j > 0 && j < my - 1)) {
189:         value = .5 * two * (hydhx + hxdhy);
190:         PetscCall(MatSetValues(jac, 1, &grow, 1, &grow, &value, INSERT_VALUES));
191:       } else {
192:         value = .25 * two * (hydhx + hxdhy);
193:         PetscCall(MatSetValues(jac, 1, &grow, 1, &grow, &value, INSERT_VALUES));
194:       }
195:     }
196:   }
197:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &ltog));
198:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
199:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
200:   PetscFunctionReturn(PETSC_SUCCESS);
201: }

203: /*TEST

205:     test:
206:       args: -ksp_monitor_short

208:     test:
209:       suffix: 2
210:       args: -ksp_monitor_short
211:       nsize: 3

213:     test:
214:       suffix: ml_1
215:       args: -ksp_monitor_short -pc_type ml -mat_no_inode
216:       nsize: 3
217:       requires: ml

219:     test:
220:       suffix: ml_2
221:       args: -ksp_monitor_short -pc_type ml -mat_no_inode -ksp_max_it 3
222:       nsize: 3
223:       requires: ml

225:     test:
226:       suffix: ml_3
227:       args: -ksp_monitor_short -pc_type ml -mat_no_inode -pc_mg_type ADDITIVE -ksp_max_it 7
228:       nsize: 1
229:       requires: ml

231:     test:
232:       suffix: cycles
233:       nsize: {{1 2}}
234:       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

236:     test:
237:       suffix: matcycles
238:       nsize: {{1 2}}
239:       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}

241:     test:
242:       requires: ml
243:       suffix: matcycles_ml
244:       nsize: {{1 2}}
245:       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}

247:     testset:
248:       requires: hpddm
249:       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
250:       test:
251:         suffix: matcycles_hpddm_mg
252:         nsize: {{1 2}}
253:         args: -pc_mg_type {{additive multiplicative full kaskade}separate output} -ksp_matsolve_batch_size {{4 7}separate output}
254:       test:
255:         requires: !__float128 !__fp16
256:         suffix: hpddm_mg_mixed_precision
257:         nsize: 2
258:         output_file: output/ex26_matcycles_hpddm_mg_pc_mg_type-multiplicative_ksp_matsolve_batch_size-4.out
259:         args: -ksp_matsolve_batch_size 4 -ksp_hpddm_precision {{single double}shared output}
260:       test:
261:         requires: __float128
262:         suffix: hpddm_mg_mixed_precision___float128
263:         nsize: 2
264:         output_file: output/ex26_matcycles_hpddm_mg_pc_mg_type-multiplicative_ksp_matsolve_batch_size-4.out
265:         args: -ksp_matsolve_batch_size 4 -ksp_hpddm_precision {{double quadruple}shared output}
266:       test:
267:         requires: double defined(PETSC_HAVE_F2CBLASLAPACK___FLOAT128_BINDINGS)
268:         suffix: hpddm_mg_mixed_precision_double
269:         nsize: 2
270:         output_file: output/ex26_matcycles_hpddm_mg_pc_mg_type-multiplicative_ksp_matsolve_batch_size-4.out
271:         args: -ksp_matsolve_batch_size 4 -ksp_hpddm_precision quadruple
272:       test:
273:         requires: single defined(PETSC_HAVE_F2CBLASLAPACK___FP16_BINDINGS)
274:         suffix: hpddm_mg_mixed_precision_single
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 half -ksp_rtol 1e-3

279:     test:
280:       requires: hpddm
281:       nsize: {{1 2}}
282:       suffix: matcycles_hpddm_ilu
283:       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}

285: TEST*/