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