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;

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

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

 44:   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:         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:         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:         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:         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:       MatSetValues(A, 1, &II, 1, &II, &v, ADD_VALUES);
 86:     }
 87:   }

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

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

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

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

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

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

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

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

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

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

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

160:   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*/