Actual source code: ex4.c


  2: static char help[] = "Solves a linear system in parallel with KSP and HMG.\n\
  3: Input parameters include:\n\
  4:   -view_exact_sol    : write exact solution vector to stdout\n\
  5:   -m  <mesh_x>       : number of mesh points in x-direction\n\
  6:   -n  <mesh_y>       : number of mesh points in y-direction\n\
  7:   -bs                : number of variables on each mesh vertex \n\n";

  9: /*
 10:   Simple example is used to test PCHMG
 11: */
 12: #include <petscksp.h>

 14: int main(int argc, char **args)
 15: {
 16:   Vec         x, b, u; /* approx solution, RHS, exact solution */
 17:   Mat         A;       /* linear system matrix */
 18:   KSP         ksp;     /* linear solver context */
 19:   PetscReal   norm;    /* norm of solution error */
 20:   PetscInt    i, j, Ii, J, Istart, Iend, m = 8, n = 7, its, bs = 1, II, JJ, jj;
 21:   PetscBool   flg, test = PETSC_FALSE, reuse = PETSC_FALSE, viewexpl = PETSC_FALSE;
 22:   PetscScalar v;
 23:   PC          pc;

 25:   PetscFunctionBeginUser;
 26:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 27:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 28:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 29:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
 30:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_hmg_interface", &test, NULL));
 31:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_reuse_interpolation", &reuse, NULL));
 32:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_explicit_mat", &viewexpl, NULL));

 34:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 35:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n * bs, m * n * bs));
 36:   PetscCall(MatSetBlockSize(A, bs));
 37:   PetscCall(MatSetFromOptions(A));
 38:   PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
 39:   PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
 40: #if defined(PETSC_HAVE_HYPRE)
 41:   PetscCall(MatHYPRESetPreallocation(A, 5, NULL, 5, NULL));
 42: #endif

 44:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

 46:   for (Ii = Istart / bs; Ii < Iend / bs; Ii++) {
 47:     v = -1.0;
 48:     i = Ii / n;
 49:     j = Ii - i * n;
 50:     if (i > 0) {
 51:       J = Ii - n;
 52:       for (jj = 0; jj < bs; jj++) {
 53:         II = Ii * bs + jj;
 54:         JJ = J * bs + jj;
 55:         PetscCall(MatSetValues(A, 1, &II, 1, &JJ, &v, ADD_VALUES));
 56:       }
 57:     }
 58:     if (i < m - 1) {
 59:       J = Ii + n;
 60:       for (jj = 0; jj < bs; jj++) {
 61:         II = Ii * bs + jj;
 62:         JJ = J * bs + jj;
 63:         PetscCall(MatSetValues(A, 1, &II, 1, &JJ, &v, ADD_VALUES));
 64:       }
 65:     }
 66:     if (j > 0) {
 67:       J = Ii - 1;
 68:       for (jj = 0; jj < bs; jj++) {
 69:         II = Ii * bs + jj;
 70:         JJ = J * bs + jj;
 71:         PetscCall(MatSetValues(A, 1, &II, 1, &JJ, &v, ADD_VALUES));
 72:       }
 73:     }
 74:     if (j < n - 1) {
 75:       J = Ii + 1;
 76:       for (jj = 0; jj < bs; jj++) {
 77:         II = Ii * bs + jj;
 78:         JJ = J * bs + jj;
 79:         PetscCall(MatSetValues(A, 1, &II, 1, &JJ, &v, ADD_VALUES));
 80:       }
 81:     }
 82:     v = 4.0;
 83:     for (jj = 0; jj < bs; jj++) {
 84:       II = Ii * bs + jj;
 85:       PetscCall(MatSetValues(A, 1, &II, 1, &II, &v, ADD_VALUES));
 86:     }
 87:   }

 89:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 90:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 91:   if (viewexpl) {
 92:     Mat E;
 93:     PetscCall(MatComputeOperator(A, MATAIJ, &E));
 94:     PetscCall(MatView(E, NULL));
 95:     PetscCall(MatDestroy(&E));
 96:   }

 98:   PetscCall(MatCreateVecs(A, &u, NULL));
 99:   PetscCall(VecSetFromOptions(u));
100:   PetscCall(VecDuplicate(u, &b));
101:   PetscCall(VecDuplicate(b, &x));

103:   PetscCall(VecSet(u, 1.0));
104:   PetscCall(MatMult(A, u, b));

106:   flg = PETSC_FALSE;
107:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_exact_sol", &flg, NULL));
108:   if (flg) PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));

110:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
111:   PetscCall(KSPSetOperators(ksp, A, A));
112:   PetscCall(KSPSetTolerances(ksp, 1.e-2 / ((m + 1) * (n + 1)), 1.e-50, PETSC_DEFAULT, PETSC_DEFAULT));

114:   if (test) {
115:     PetscCall(KSPGetPC(ksp, &pc));
116:     PetscCall(PCSetType(pc, PCHMG));
117:     PetscCall(PCHMGSetInnerPCType(pc, PCGAMG));
118:     PetscCall(PCHMGSetReuseInterpolation(pc, PETSC_TRUE));
119:     PetscCall(PCHMGSetUseSubspaceCoarsening(pc, PETSC_TRUE));
120:     PetscCall(PCHMGUseMatMAIJ(pc, PETSC_FALSE));
121:     PetscCall(PCHMGSetCoarseningComponent(pc, 0));
122:   }

124:   PetscCall(KSPSetFromOptions(ksp));
125:   PetscCall(KSPSolve(ksp, b, x));

127:   if (reuse) {
128:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
129:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
130:     PetscCall(KSPSolve(ksp, b, x));
131:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
132:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
133:     PetscCall(KSPSolve(ksp, b, x));
134:     /* Make sparsity pattern different and reuse interpolation */
135:     PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
136:     PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
137:     PetscCall(MatGetSize(A, &m, NULL));
138:     n = 0;
139:     v = 0;
140:     m--;
141:     /* Connect the last element to the first element */
142:     PetscCall(MatSetValue(A, m, n, v, ADD_VALUES));
143:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
144:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
145:     PetscCall(KSPSolve(ksp, b, x));
146:   }

148:   PetscCall(VecAXPY(x, -1.0, u));
149:   PetscCall(VecNorm(x, NORM_2, &norm));
150:   PetscCall(KSPGetIterationNumber(ksp, &its));

152:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));

154:   PetscCall(KSPDestroy(&ksp));
155:   PetscCall(VecDestroy(&u));
156:   PetscCall(VecDestroy(&x));
157:   PetscCall(VecDestroy(&b));
158:   PetscCall(MatDestroy(&A));

160:   PetscCall(PetscFinalize());
161:   return 0;
162: }

164: /*TEST

166:    build:
167:       requires: !complex !single

169:    test:
170:       suffix: hypre
171:       nsize: 2
172:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
173:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre

175:    test:
176:       suffix: hypre_bs4
177:       nsize: 2
178:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
179:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1

181:    test:
182:       suffix: hypre_asm
183:       nsize: 2
184:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
185:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1 -mg_levels_3_pc_type asm

187:    test:
188:       suffix: hypre_fieldsplit
189:       nsize: 2
190:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
191:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -mg_levels_4_pc_type fieldsplit

193:    test:
194:       suffix: gamg
195:       nsize: 2
196:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg

198:    test:
199:       suffix: gamg_bs4
200:       nsize: 2
201:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1

203:    test:
204:       suffix: gamg_asm
205:       nsize: 2
206:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1 -mg_levels_1_pc_type asm

208:    test:
209:       suffix: gamg_fieldsplit
210:       nsize: 2
211:       args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -mg_levels_1_pc_type fieldsplit

213:    test:
214:       suffix: interface
215:       nsize: 2
216:       args: -ksp_monitor -ksp_rtol 1e-6 -test_hmg_interface 1 -bs 4

218:    test:
219:       suffix: reuse
220:       nsize: 2
221:       args: -ksp_monitor -ksp_rtol 1e-6   -pc_type hmg -pc_hmg_reuse_interpolation 1 -test_reuse_interpolation 1 -hmg_inner_pc_type gamg

223:    test:
224:       suffix: component
225:       nsize: 2
226:       args: -ksp_monitor -ksp_rtol 1e-6 -pc_type hmg -pc_hmg_coarsening_component 2  -pc_hmg_use_subspace_coarsening 1 -bs 4 -hmg_inner_pc_type gamg

228:    testset:
229:       output_file: output/ex4_expl.out
230:       nsize: {{1 2}}
231:       filter: grep -v " MPI process" | grep -v " type:" | grep -v "Mat Object"
232:       args: -ksp_converged_reason -view_explicit_mat -pc_type none -ksp_type {{cg gmres}}
233:       test:
234:         suffix: expl_aij
235:         args: -mat_type aij
236:       test:
237:         suffix: expl_hypre
238:         requires: hypre
239:         args: -mat_type hypre

241:    test:
242:       suffix: hypre_device
243:       nsize: {{1 2}}
244:       requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
245:       args: -mat_type hypre -ksp_converged_reason -pc_type hypre -m 13 -n 17

247:    test:
248:       suffix: hypre_device_cusparse
249:       output_file: output/ex4_hypre_device.out
250:       nsize: {{1 2}}
251:       requires: hypre cuda defined(PETSC_HAVE_HYPRE_DEVICE)
252:       args: -mat_type {{aij aijcusparse}} -vec_type cuda -ksp_converged_reason -pc_type hypre -m 13 -n 17

254: TEST*/