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 x,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;
23: PetscErrorCode ierr;
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,"-pc_hpddm_block_splitting",&flg,NULL);
82: if (!flg) {
83: PCHPDDMSetAuxiliaryMat(pc,is,aux,NULL,NULL);
84: PCHPDDMHasNeumannMat(pc,PETSC_FALSE); /* PETSC_TRUE is fine as well, just testing */
85: }
86: flg = PETSC_FALSE;
87: PetscOptionsGetBool(NULL,NULL,"-set_rhs",&flg,NULL);
88: if (flg) { /* user-provided RHS for concurrent generalized eigenvalue problems */
89: Mat a,c,P; /* usually assembled automatically in PCHPDDM, this is solely for testing PCHPDDMSetRHSMat() */
90: PetscInt rstart,rend,location;
91: MatDuplicate(aux,MAT_DO_NOT_COPY_VALUES,&B); /* duplicate so that MatStructure is SAME_NONZERO_PATTERN */
92: MatGetDiagonalBlock(A,&a);
93: MatGetOwnershipRange(A,&rstart,&rend);
94: ISGetLocalSize(is,&m);
95: MatCreateSeqAIJ(PETSC_COMM_SELF,rend-rstart,m,1,NULL,&P);
96: for (m = rstart; m < rend; ++m) {
97: ISLocate(is,m,&location);
99: MatSetValue(P,m-rstart,location,1.0,INSERT_VALUES);
100: }
101: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
102: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
103: PetscObjectTypeCompare((PetscObject)a,MATSEQAIJ,&flg);
104: if (flg) {
105: MatPtAP(a,P,MAT_INITIAL_MATRIX,1.0,&X); /* MatPtAP() is used to extend diagonal blocks with zeros on the overlap */
106: } else { /* workaround for MatPtAP() limitations with some types */
107: MatConvert(a,MATSEQAIJ,MAT_INITIAL_MATRIX,&c);
108: MatPtAP(c,P,MAT_INITIAL_MATRIX,1.0,&X);
109: MatDestroy(&c);
110: }
111: MatDestroy(&P);
112: MatAXPY(B,1.0,X,SUBSET_NONZERO_PATTERN);
113: MatDestroy(&X);
114: MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
115: PCHPDDMSetRHSMat(pc,B);
116: MatDestroy(&B);
117: }
118: #endif
119: ISDestroy(&is);
120: MatDestroy(&aux);
121: KSPSetFromOptions(ksp);
122: MatCreateVecs(A,&x,&b);
123: VecSet(b,1.0);
124: KSPSolve(ksp,b,x);
125: VecGetLocalSize(x,&m);
126: VecDestroy(&x);
127: VecDestroy(&b);
128: if (N > 1) {
129: KSPType type;
130: PetscOptionsClearValue(NULL,"-ksp_converged_reason");
131: KSPSetFromOptions(ksp);
132: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&B);
133: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&X);
134: MatSetRandom(B,NULL);
135: /* this is algorithmically optimal in the sense that blocks of vectors are coarsened or interpolated using matrix--matrix operations */
136: /* PCHPDDM however heavily relies on MPI[S]BAIJ format for which there is no efficient MatProduct implementation */
137: KSPMatSolve(ksp,B,X);
138: KSPGetType(ksp,&type);
139: PetscStrcmp(type,KSPHPDDM,&flg);
140: #if defined(PETSC_HAVE_HPDDM)
141: if (flg) {
142: PetscReal norm;
143: KSPHPDDMType type;
144: KSPHPDDMGetType(ksp,&type);
145: if (type == KSP_HPDDM_TYPE_PREONLY || type == KSP_HPDDM_TYPE_CG || type == KSP_HPDDM_TYPE_GMRES || type == KSP_HPDDM_TYPE_GCRODR) {
146: Mat C;
147: MatDuplicate(X,MAT_DO_NOT_COPY_VALUES,&C);
148: KSPSetMatSolveBatchSize(ksp,1);
149: KSPMatSolve(ksp,B,C);
150: MatAYPX(C,-1.0,X,SAME_NONZERO_PATTERN);
151: MatNorm(C,NORM_INFINITY,&norm);
152: MatDestroy(&C);
154: }
155: }
156: #endif
157: MatDestroy(&X);
158: MatDestroy(&B);
159: }
160: PetscObjectTypeCompare((PetscObject)pc,PCHPDDM,&flg);
161: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
162: if (flg) {
163: PCHPDDMGetSTShareSubKSP(pc,&flg);
164: }
165: #endif
166: if (flg && PetscDefined(USE_LOG)) {
167: PetscLogEventRegister("MatLUFactorSym",PC_CLASSID,&event);
168: PetscLogEventGetPerfInfo(PETSC_DETERMINE,event,&info1);
169: PetscLogEventRegister("MatLUFactorNum",PC_CLASSID,&event);
170: PetscLogEventGetPerfInfo(PETSC_DETERMINE,event,&info2);
171: if (info1.count || info2.count) {
173: } else {
174: PetscLogEventRegister("MatCholFctrSym",PC_CLASSID,&event);
175: PetscLogEventGetPerfInfo(PETSC_DETERMINE,event,&info1);
176: PetscLogEventRegister("MatCholFctrNum",PC_CLASSID,&event);
177: PetscLogEventGetPerfInfo(PETSC_DETERMINE,event,&info2);
179: }
180: }
181: KSPDestroy(&ksp);
182: MatDestroy(&A);
183: PetscFinalize();
184: return 0;
185: }
187: /*TEST
189: test:
190: requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
191: nsize: 4
192: 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
194: testset:
195: requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
196: nsize: 4
197: 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
198: test:
199: suffix: geneo
200: 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}
201: test:
202: suffix: geneo_block_splitting
203: output_file: output/ex76_geneo_pc_hpddm_levels_1_eps_nev-15.out
204: 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"
205: 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}
206: test:
207: suffix: geneo_share
208: output_file: output/ex76_geneo_pc_hpddm_levels_1_eps_nev-5.out
209: args: -pc_hpddm_levels_1_st_pc_type cholesky -pc_hpddm_levels_1_eps_nev 5 -pc_hpddm_levels_1_st_share_sub_ksp
211: testset:
212: requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
213: nsize: 4
214: 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
215: test:
216: suffix: geneo_share_cholesky
217: output_file: output/ex76_geneo_share.out
218: # extra -pc_hpddm_levels_1_eps_gen_non_hermitian needed to avoid failures with PETSc Cholesky
219: 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}
220: test:
221: suffix: geneo_share_cholesky_matstructure
222: output_file: output/ex76_geneo_share.out
223: # extra -pc_hpddm_levels_1_eps_gen_non_hermitian needed to avoid failures with PETSc Cholesky
224: 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}
225: test:
226: requires: mumps
227: suffix: geneo_share_lu
228: output_file: output/ex76_geneo_share.out
229: # extra -pc_factor_mat_solver_type mumps needed to avoid failures with PETSc LU
230: 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}
231: test:
232: requires: mumps
233: suffix: geneo_share_lu_matstructure
234: output_file: output/ex76_geneo_share.out
235: # extra -pc_factor_mat_solver_type mumps needed to avoid failures with PETSc LU
236: 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
238: test:
239: requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
240: suffix: fgmres_geneo_20_p_2
241: nsize: 4
242: 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
244: testset:
245: requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
246: output_file: output/ex76_fgmres_geneo_20_p_2.out
247: nsize: 4
248: 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
249: test:
250: suffix: fgmres_geneo_20_p_2_geneo
251: args: -mat_type {{aij sbaij}shared output}
252: test:
253: suffix: fgmres_geneo_20_p_2_geneo_algebraic
254: args: -pc_hpddm_levels_2_st_pc_type mat
255: # PCHPDDM + KSPHPDDM test to exercise multilevel + multiple RHS in one go
256: test:
257: requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
258: suffix: fgmres_geneo_20_p_2_geneo_rhs
259: output_file: output/ex76_fgmres_geneo_20_p_2.out
260: # for -pc_hpddm_coarse_correction additive
261: filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 37/Linear solve converged due to CONVERGED_RTOL iterations 25/g"
262: nsize: 4
263: 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}
265: testset:
266: 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)
267: filter: egrep -e "Linear solve" -e " executing" | sed -e "s/MPI = 1/MPI = 2/g" -e "s/OMP = 1/OMP = 2/g"
268: nsize: 4
269: 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}
270: test:
271: suffix: geneo_mumps_use_omp_threads_1
272: output_file: output/ex76_geneo_mumps_use_omp_threads.out
273: args: -pc_hpddm_coarse_mat_type {{baij sbaij}shared output}
274: test:
275: suffix: geneo_mumps_use_omp_threads_2
276: output_file: output/ex76_geneo_mumps_use_omp_threads.out
277: args: -pc_hpddm_coarse_mat_type aij -pc_hpddm_levels_1_eps_threshold 0.3 -pc_hpddm_coarse_pc_type cholesky
279: test:
280: requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
281: suffix: reuse_symbolic
282: output_file: output/ex77_preonly.out
283: nsize: 4
284: 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
286: TEST*/