Actual source code: ex76.c

  1: #include <petscksp.h>

  3: static char help[] = "Solves a linear system using PCHPDDM.\n\n";

  5: int main(int argc, char **args)
  6: {
  7:   Vec             b;            /* computed solution and RHS */
  8:   Mat             A, aux, X, B; /* linear system matrix */
  9:   KSP             ksp;          /* linear solver context */
 10:   PC              pc;
 11:   IS              is, sizes;
 12:   const PetscInt *idx;
 13:   PetscMPIInt     rank, size;
 14:   PetscInt        m, N = 1;
 15:   const char     *deft = MATAIJ;
 16:   PetscViewer     viewer;
 17:   char            dir[PETSC_MAX_PATH_LEN], name[PETSC_MAX_PATH_LEN], type[256];
 18:   PetscBool       flg;
 19: #if defined(PETSC_USE_LOG)
 20:   PetscLogEvent event;
 21: #endif
 22:   PetscEventPerfInfo info1, info2;

 25:   PetscInitialize(&argc, &args, NULL, help);
 26:   PetscLogDefaultBegin();
 27:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 29:   PetscOptionsGetInt(NULL, NULL, "-rhs", &N, NULL);
 30:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 31:   MatCreate(PETSC_COMM_WORLD, &A);
 32:   MatCreate(PETSC_COMM_SELF, &aux);
 33:   ISCreate(PETSC_COMM_SELF, &is);
 34:   PetscStrcpy(dir, ".");
 35:   PetscOptionsGetString(NULL, NULL, "-load_dir", dir, sizeof(dir), NULL);
 36:   /* loading matrices */
 37:   PetscSNPrintf(name, sizeof(name), "%s/sizes_%d_%d.dat", dir, rank, size);
 38:   PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_READ, &viewer);
 39:   ISCreate(PETSC_COMM_SELF, &sizes);
 40:   ISLoad(sizes, viewer);
 41:   ISGetIndices(sizes, &idx);
 42:   MatSetSizes(A, idx[0], idx[1], idx[2], idx[3]);
 43:   ISRestoreIndices(sizes, &idx);
 44:   ISDestroy(&sizes);
 45:   PetscViewerDestroy(&viewer);
 46:   MatSetUp(A);
 47:   PetscSNPrintf(name, sizeof(name), "%s/A.dat", dir);
 48:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer);
 49:   MatLoad(A, viewer);
 50:   PetscViewerDestroy(&viewer);
 51:   PetscSNPrintf(name, sizeof(name), "%s/is_%d_%d.dat", dir, rank, size);
 52:   PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_READ, &viewer);
 53:   ISLoad(is, viewer);
 54:   PetscViewerDestroy(&viewer);
 55:   PetscSNPrintf(name, sizeof(name), "%s/Neumann_%d_%d.dat", dir, rank, size);
 56:   PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_READ, &viewer);
 57:   MatLoad(aux, viewer);
 58:   PetscViewerDestroy(&viewer);
 59:   flg = PETSC_FALSE;
 60:   PetscOptionsGetBool(NULL, NULL, "-pc_hpddm_levels_1_st_share_sub_ksp", &flg, NULL);
 61:   if (flg) { /* PETSc LU/Cholesky is struggling numerically for bs > 1          */
 62:              /* only set the proper bs for the geneo_share_* tests, 1 otherwise */
 63:     MatSetBlockSizesFromMats(aux, A, A);
 64:   }
 65:   MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
 66:   MatSetOption(aux, MAT_SYMMETRIC, PETSC_TRUE);
 67:   /* ready for testing */
 68:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
 69:   PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, deft, type, 256, &flg);
 70:   PetscOptionsEnd();
 71:   if (flg) {
 72:     MatConvert(A, type, MAT_INPLACE_MATRIX, &A);
 73:     MatConvert(aux, type, MAT_INPLACE_MATRIX, &aux);
 74:   }
 75:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 76:   KSPSetOperators(ksp, A, A);
 77:   KSPGetPC(ksp, &pc);
 78:   PCSetType(pc, PCHPDDM);
 79: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
 80:   flg = PETSC_FALSE;
 81:   PetscOptionsGetBool(NULL, NULL, "-reset", &flg, NULL);
 82:   if (flg) {
 83:     PetscOptionsSetValue(NULL, "-pc_hpddm_block_splitting", "true");
 84:     PCSetFromOptions(pc);
 85:     PCSetUp(pc);
 86:     PetscOptionsClearValue(NULL, "-pc_hpddm_block_splitting");
 87:   }
 88:   PCHPDDMSetAuxiliaryMat(pc, is, aux, NULL, NULL);
 89:   PCHPDDMHasNeumannMat(pc, PETSC_FALSE); /* PETSC_TRUE is fine as well, just testing */
 90:   flg = PETSC_FALSE;
 91:   PetscOptionsGetBool(NULL, NULL, "-set_rhs", &flg, NULL);
 92:   if (flg) {          /* user-provided RHS for concurrent generalized eigenvalue problems                                 */
 93:     Mat      a, c, P; /* usually assembled automatically in PCHPDDM, this is solely for testing PCHPDDMSetRHSMat() */
 94:     PetscInt rstart, rend, location;
 95:     MatDuplicate(aux, MAT_DO_NOT_COPY_VALUES, &B); /* duplicate so that MatStructure is SAME_NONZERO_PATTERN */
 96:     MatGetDiagonalBlock(A, &a);
 97:     MatGetOwnershipRange(A, &rstart, &rend);
 98:     ISGetLocalSize(is, &m);
 99:     MatCreateSeqAIJ(PETSC_COMM_SELF, rend - rstart, m, 1, NULL, &P);
100:     for (m = rstart; m < rend; ++m) {
101:       ISLocate(is, m, &location);
103:       MatSetValue(P, m - rstart, location, 1.0, INSERT_VALUES);
104:     }
105:     MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
106:     MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
107:     PetscObjectTypeCompare((PetscObject)a, MATSEQAIJ, &flg);
108:     if (flg) MatPtAP(a, P, MAT_INITIAL_MATRIX, 1.0, &X); /* MatPtAP() is used to extend diagonal blocks with zeros on the overlap */
109:     else { /* workaround for MatPtAP() limitations with some types */ MatConvert(a, MATSEQAIJ, MAT_INITIAL_MATRIX, &c);
110:       MatPtAP(c, P, MAT_INITIAL_MATRIX, 1.0, &X);
111:       MatDestroy(&c);
112:     }
113:     MatDestroy(&P);
114:     MatAXPY(B, 1.0, X, SUBSET_NONZERO_PATTERN);
115:     MatDestroy(&X);
116:     MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE);
117:     PCHPDDMSetRHSMat(pc, B);
118:     MatDestroy(&B);
119:   }
120: #endif
121:   MatDestroy(&aux);
122:   KSPSetFromOptions(ksp);
123:   PetscObjectTypeCompare((PetscObject)pc, PCASM, &flg);
124:   if (flg) {
125:     flg = PETSC_FALSE;
126:     PetscOptionsGetBool(NULL, NULL, "-pc_hpddm_define_subdomains", &flg, NULL);
127:     if (flg) {
128:       IS rows;
129:       MatGetOwnershipIS(A, &rows, NULL);
130:       PCASMSetLocalSubdomains(pc, 1, &is, &rows);
131:       ISDestroy(&rows);
132:     }
133:   }
134:   ISDestroy(&is);
135:   MatCreateVecs(A, NULL, &b);
136:   VecSet(b, 1.0);
137:   KSPSolve(ksp, b, b);
138:   VecGetLocalSize(b, &m);
139:   VecDestroy(&b);
140:   if (N > 1) {
141:     KSPType type;
142:     PetscOptionsClearValue(NULL, "-ksp_converged_reason");
143:     KSPSetFromOptions(ksp);
144:     MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B);
145:     MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &X);
146:     MatSetRandom(B, NULL);
147:     /* this is algorithmically optimal in the sense that blocks of vectors are coarsened or interpolated using matrix--matrix operations */
148:     /* PCHPDDM however heavily relies on MPI[S]BAIJ format for which there is no efficient MatProduct implementation */
149:     KSPMatSolve(ksp, B, X);
150:     KSPGetType(ksp, &type);
151:     PetscStrcmp(type, KSPHPDDM, &flg);
152: #if defined(PETSC_HAVE_HPDDM)
153:     if (flg) {
154:       PetscReal    norm;
155:       KSPHPDDMType type;
156:       KSPHPDDMGetType(ksp, &type);
157:       if (type == KSP_HPDDM_TYPE_PREONLY || type == KSP_HPDDM_TYPE_CG || type == KSP_HPDDM_TYPE_GMRES || type == KSP_HPDDM_TYPE_GCRODR) {
158:         Mat C;
159:         MatDuplicate(X, MAT_DO_NOT_COPY_VALUES, &C);
160:         KSPSetMatSolveBatchSize(ksp, 1);
161:         KSPMatSolve(ksp, B, C);
162:         MatAYPX(C, -1.0, X, SAME_NONZERO_PATTERN);
163:         MatNorm(C, NORM_INFINITY, &norm);
164:         MatDestroy(&C);
166:       }
167:     }
168: #endif
169:     MatDestroy(&X);
170:     MatDestroy(&B);
171:   }
172:   PetscObjectTypeCompare((PetscObject)pc, PCHPDDM, &flg);
173: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
174:   if (flg) PCHPDDMGetSTShareSubKSP(pc, &flg);
175: #endif
176:   if (flg && PetscDefined(USE_LOG)) {
177:     PetscLogEventRegister("MatLUFactorSym", PC_CLASSID, &event);
178:     PetscLogEventGetPerfInfo(PETSC_DETERMINE, event, &info1);
179:     PetscLogEventRegister("MatLUFactorNum", PC_CLASSID, &event);
180:     PetscLogEventGetPerfInfo(PETSC_DETERMINE, event, &info2);
181:     if (!info1.count && !info2.count) {
182:       PetscLogEventRegister("MatCholFctrSym", PC_CLASSID, &event);
183:       PetscLogEventGetPerfInfo(PETSC_DETERMINE, event, &info1);
184:       PetscLogEventRegister("MatCholFctrNum", PC_CLASSID, &event);
185:       PetscLogEventGetPerfInfo(PETSC_DETERMINE, event, &info2);
188:   }
189:   KSPDestroy(&ksp);
190:   MatDestroy(&A);
191:   PetscFinalize();
192:   return 0;
193: }

195: /*TEST

197:    test:
198:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
199:       nsize: 4
200:       args: -ksp_rtol 1e-3 -ksp_converged_reason -pc_type {{bjacobi hpddm}shared output} -pc_hpddm_coarse_sub_pc_type lu -sub_pc_type lu -options_left no -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO

202:    test:
203:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
204:       suffix: define_subdomains
205:       nsize: 4
206:       args: -ksp_rtol 1e-3 -ksp_converged_reason -pc_type {{asm hpddm}shared output} -pc_hpddm_coarse_sub_pc_type lu -sub_pc_type lu -pc_hpddm_define_subdomains -options_left no -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO

208:    testset:
209:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
210:       nsize: 4
211:       args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_coarse_pc_type redundant -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
212:       test:
213:         suffix: geneo
214:         args: -pc_hpddm_coarse_p {{1 2}shared output} -pc_hpddm_levels_1_st_pc_type cholesky -pc_hpddm_levels_1_eps_nev {{5 15}separate output} -mat_type {{aij baij sbaij}shared output}
215:       test:
216:         suffix: geneo_block_splitting
217:         output_file: output/ex76_geneo_pc_hpddm_levels_1_eps_nev-15.out
218:         filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 1[6-9]/Linear solve converged due to CONVERGED_RTOL iterations 11/g"
219:         args: -pc_hpddm_coarse_p 2 -pc_hpddm_levels_1_eps_nev 15 -pc_hpddm_block_splitting -pc_hpddm_levels_1_st_pc_type lu -pc_hpddm_levels_1_eps_gen_non_hermitian -mat_type {{aij baij}shared output}
220:       test:
221:         suffix: geneo_share
222:         output_file: output/ex76_geneo_pc_hpddm_levels_1_eps_nev-5.out
223:         args: -pc_hpddm_levels_1_st_pc_type cholesky -pc_hpddm_levels_1_eps_nev 5 -pc_hpddm_levels_1_st_share_sub_ksp -reset {{false true}shared output}

225:    testset:
226:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
227:       nsize: 4
228:       args: -ksp_converged_reason -ksp_max_it 150 -pc_type hpddm -pc_hpddm_levels_1_eps_nev 5 -pc_hpddm_coarse_p 1 -pc_hpddm_coarse_pc_type redundant -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -pc_hpddm_define_subdomains
229:       test:
230:         suffix: geneo_share_cholesky
231:         output_file: output/ex76_geneo_share.out
232:         # extra -pc_hpddm_levels_1_eps_gen_non_hermitian needed to avoid failures with PETSc Cholesky
233:         args: -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_st_pc_type cholesky -mat_type {{aij baij sbaij}shared output} -pc_hpddm_levels_1_eps_gen_non_hermitian -pc_hpddm_has_neumann -pc_hpddm_levels_1_st_share_sub_ksp {{false true}shared output}
234:       test:
235:         suffix: geneo_share_cholesky_matstructure
236:         output_file: output/ex76_geneo_share.out
237:         # extra -pc_hpddm_levels_1_eps_gen_non_hermitian needed to avoid failures with PETSc Cholesky
238:         args: -pc_hpddm_levels_1_sub_pc_type cholesky -mat_type {{baij sbaij}shared output} -pc_hpddm_levels_1_eps_gen_non_hermitian -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_st_matstructure same -set_rhs {{false true} shared output}
239:       test:
240:         requires: mumps
241:         suffix: geneo_share_lu
242:         output_file: output/ex76_geneo_share.out
243:         # extra -pc_factor_mat_solver_type mumps needed to avoid failures with PETSc LU
244:         args: -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_st_pc_type lu -mat_type baij -pc_hpddm_levels_1_st_pc_factor_mat_solver_type mumps -pc_hpddm_levels_1_sub_pc_factor_mat_solver_type mumps -pc_hpddm_has_neumann -pc_hpddm_levels_1_st_share_sub_ksp {{false true}shared output}
245:       test:
246:         requires: mumps
247:         suffix: geneo_share_lu_matstructure
248:         output_file: output/ex76_geneo_share.out
249:         # extra -pc_factor_mat_solver_type mumps needed to avoid failures with PETSc LU
250:         args: -pc_hpddm_levels_1_sub_pc_type lu -mat_type baij -pc_hpddm_levels_1_sub_pc_factor_mat_solver_type mumps -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_st_matstructure {{same different}shared output} -pc_hpddm_levels_1_st_pc_type lu -pc_hpddm_levels_1_st_pc_factor_mat_solver_type mumps
251:       test:
252:         suffix: geneo_share_not_asm
253:         output_file: output/ex76_geneo_pc_hpddm_levels_1_eps_nev-5.out
254:         # extra -pc_hpddm_levels_1_eps_gen_non_hermitian needed to avoid failures with PETSc Cholesky
255:         args: -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_st_pc_type cholesky -pc_hpddm_levels_1_eps_gen_non_hermitian -pc_hpddm_has_neumann -pc_hpddm_levels_1_st_share_sub_ksp true -pc_hpddm_levels_1_pc_type gasm

257:    test:
258:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
259:       suffix: fgmres_geneo_20_p_2
260:       nsize: 4
261:       args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_pc_type redundant -ksp_type fgmres -pc_hpddm_coarse_mat_type {{baij sbaij}shared output} -pc_hpddm_log_separate {{false true}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO

263:    testset:
264:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
265:       output_file: output/ex76_fgmres_geneo_20_p_2.out
266:       nsize: 4
267:       args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_2_p 2 -pc_hpddm_levels_2_mat_type {{baij sbaij}shared output} -pc_hpddm_levels_2_eps_nev {{5 20}shared output} -pc_hpddm_levels_2_sub_pc_type cholesky -pc_hpddm_levels_2_ksp_type gmres -ksp_type fgmres -pc_hpddm_coarse_mat_type {{baij sbaij}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
268:       test:
269:         suffix: fgmres_geneo_20_p_2_geneo
270:         args: -mat_type {{aij sbaij}shared output}
271:       test:
272:         suffix: fgmres_geneo_20_p_2_geneo_algebraic
273:         args: -pc_hpddm_levels_2_st_pc_type mat
274:    # PCHPDDM + KSPHPDDM test to exercise multilevel + multiple RHS in one go
275:    test:
276:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
277:       suffix: fgmres_geneo_20_p_2_geneo_rhs
278:       output_file: output/ex76_fgmres_geneo_20_p_2.out
279:       # for -pc_hpddm_coarse_correction additive
280:       filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 37/Linear solve converged due to CONVERGED_RTOL iterations 25/g"
281:       nsize: 4
282:       args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_2_p 2 -pc_hpddm_levels_2_mat_type baij -pc_hpddm_levels_2_eps_nev 5 -pc_hpddm_levels_2_sub_pc_type cholesky -pc_hpddm_levels_2_ksp_max_it 10 -pc_hpddm_levels_2_ksp_type hpddm -pc_hpddm_levels_2_ksp_hpddm_type gmres -ksp_type hpddm -ksp_hpddm_variant flexible -pc_hpddm_coarse_mat_type baij -mat_type aij -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -rhs 4 -pc_hpddm_coarse_correction {{additive deflated balanced}shared output}

284:    testset:
285:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES) mumps defined(PETSC_HAVE_OPENMP_SUPPORT)
286:       filter: grep -E -e "Linear solve" -e "      executing" | sed -e "s/MPI =      1/MPI =      2/g" -e "s/OMP =      1/OMP =      2/g"
287:       nsize: 4
288:       args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 15 -pc_hpddm_levels_1_st_pc_type cholesky -pc_hpddm_coarse_p {{1 2}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -pc_hpddm_coarse_pc_factor_mat_solver_type mumps -pc_hpddm_coarse_mat_mumps_icntl_4 2 -pc_hpddm_coarse_mat_mumps_use_omp_threads {{1 2}shared output}
289:       test:
290:         suffix: geneo_mumps_use_omp_threads_1
291:         output_file: output/ex76_geneo_mumps_use_omp_threads.out
292:         args: -pc_hpddm_coarse_mat_type {{baij sbaij}shared output}
293:       test:
294:         suffix: geneo_mumps_use_omp_threads_2
295:         output_file: output/ex76_geneo_mumps_use_omp_threads.out
296:         args: -pc_hpddm_coarse_mat_type aij -pc_hpddm_levels_1_eps_threshold 0.3 -pc_hpddm_coarse_pc_type cholesky

298:    test:
299:       requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
300:       suffix: reuse_symbolic
301:       output_file: output/ex77_preonly.out
302:       nsize: 4
303:       args: -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 20 -rhs 4 -pc_hpddm_coarse_correction {{additive deflated balanced}shared output} -ksp_pc_side {{left right}shared output} -ksp_max_it 20 -ksp_type hpddm -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -pc_hpddm_define_subdomains -ksp_error_if_not_converged

305: TEST*/