Actual source code: ex4.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

163: /*TEST

165:    build:
166:       requires: !complex !single

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

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

180:    test:
181:       suffix: hypre_asm
182:       nsize: 2
183:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
184:       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

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

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

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

202:    test:
203:       suffix: gamg_asm
204:       nsize: 2
205:       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

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

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

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

222:    test:
223:       suffix: component
224:       nsize: 2
225:       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

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

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

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

253: TEST*/