Actual source code: ex76.c
petsc-3.15.0 2021-04-05
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);if (ierr) return ierr;
26: PetscLogDefaultBegin();
27: MPI_Comm_size(PETSC_COMM_WORLD,&size);
28: if (size != 4) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"This example requires 4 processes");
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: PCHPDDMSetAuxiliaryMat(pc,is,aux,NULL,NULL);
81: PCHPDDMHasNeumannMat(pc,PETSC_FALSE); /* PETSC_TRUE is fine as well, just testing */
82: #endif
83: ISDestroy(&is);
84: MatDestroy(&aux);
85: KSPSetFromOptions(ksp);
86: MatCreateVecs(A,&x,&b);
87: VecSet(b,1.0);
88: KSPSolve(ksp,b,x);
89: VecGetLocalSize(x,&m);
90: VecDestroy(&x);
91: VecDestroy(&b);
92: if (N > 1) {
93: KSPType type;
94: PetscOptionsClearValue(NULL,"-ksp_converged_reason");
95: KSPSetFromOptions(ksp);
96: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&B);
97: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&X);
98: MatSetRandom(B,NULL);
99: /* this is algorithmically optimal in the sense that blocks of vectors are coarsened or interpolated using matrix--matrix operations */
100: /* PCHPDDM however heavily relies on MPI[S]BAIJ format for which there is no efficient MatProduct implementation */
101: KSPMatSolve(ksp,B,X);
102: KSPGetType(ksp,&type);
103: PetscStrcmp(type,KSPHPDDM,&flg);
104: #if defined(PETSC_HAVE_HPDDM)
105: if (flg) {
106: PetscReal norm;
107: KSPHPDDMType type;
108: KSPHPDDMGetType(ksp,&type);
109: if (type == KSP_HPDDM_TYPE_PREONLY || type == KSP_HPDDM_TYPE_CG || type == KSP_HPDDM_TYPE_GMRES || type == KSP_HPDDM_TYPE_GCRODR) {
110: Mat C;
111: MatDuplicate(X,MAT_DO_NOT_COPY_VALUES,&C);
112: KSPSetMatSolveBatchSize(ksp,1);
113: KSPMatSolve(ksp,B,C);
114: MatAYPX(C,-1.0,X,SAME_NONZERO_PATTERN);
115: MatNorm(C,NORM_INFINITY,&norm);
116: MatDestroy(&C);
117: if (norm > 100*PETSC_MACHINE_EPSILON) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"KSPMatSolve() and KSPSolve() difference has nonzero norm %g with pseudo-block KSPHPDDMType %s",(double)norm,KSPHPDDMTypes[type]);
118: }
119: }
120: #endif
121: MatDestroy(&X);
122: MatDestroy(&B);
123: }
124: PetscObjectTypeCompare((PetscObject)pc,PCHPDDM,&flg);
125: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
126: if (flg) {
127: PCHPDDMGetSTShareSubKSP(pc,&flg);
128: }
129: #endif
130: if (flg && PetscDefined(USE_LOG)) {
131: PetscLogEventRegister("MatLUFactorSym",PC_CLASSID,&event);
132: PetscLogEventGetPerfInfo(PETSC_DETERMINE,event,&info1);
133: PetscLogEventRegister("MatLUFactorNum",PC_CLASSID,&event);
134: PetscLogEventGetPerfInfo(PETSC_DETERMINE,event,&info2);
135: if (info1.count || info2.count) {
136: if (info2.count <= info1.count) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"LU numerical factorization (%d) not called more times than LU symbolic factorization (%d), broken -pc_hpddm_levels_1_st_share_sub_ksp",info2.count,info1.count);
137: } else {
138: PetscLogEventRegister("MatCholFctrSym",PC_CLASSID,&event);
139: PetscLogEventGetPerfInfo(PETSC_DETERMINE,event,&info1);
140: PetscLogEventRegister("MatCholFctrNum",PC_CLASSID,&event);
141: PetscLogEventGetPerfInfo(PETSC_DETERMINE,event,&info2);
142: if (info2.count <= info1.count) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cholesky numerical factorization (%d) not called more times than Cholesky symbolic factorization (%d), broken -pc_hpddm_levels_1_st_share_sub_ksp",info2.count,info1.count);
143: }
144: }
145: KSPDestroy(&ksp);
146: MatDestroy(&A);
147: PetscFinalize();
148: return ierr;
149: }
151: /*TEST
153: test:
154: requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_DYNAMIC_LIBRARIES) define(PETSC_USE_SHARED_LIBRARIES)
155: nsize: 4
156: 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
158: test:
159: requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_DYNAMIC_LIBRARIES) define(PETSC_USE_SHARED_LIBRARIES)
160: suffix: geneo
161: nsize: 4
162: args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev {{5 15}separate output} -pc_hpddm_levels_1_st_pc_type cholesky -pc_hpddm_coarse_p {{1 2}shared output} -pc_hpddm_coarse_pc_type redundant -mat_type {{aij baij sbaij}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
164: testset:
165: requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_DYNAMIC_LIBRARIES) define(PETSC_USE_SHARED_LIBRARIES)
166: nsize: 4
167: 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 -pc_hpddm_has_neumann -pc_hpddm_levels_1_st_share_sub_ksp {{false true}shared output}
168: test:
169: suffix: geneo_share_cholesky
170: output_file: output/ex76_geneo_share.out
171: # extra -pc_hpddm_levels_1_eps_gen_non_hermitian needed to avoid failures with PETSc Cholesky
172: 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
173: test:
174: requires: mumps
175: suffix: geneo_share_lu
176: output_file: output/ex76_geneo_share.out
177: # extra -pc_factor_mat_solver_type mumps needed to avoid failures with PETSc LU
178: 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
180: test:
181: requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_DYNAMIC_LIBRARIES) define(PETSC_USE_SHARED_LIBRARIES)
182: suffix: fgmres_geneo_20_p_2
183: nsize: 4
184: 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
186: test:
187: requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_DYNAMIC_LIBRARIES) define(PETSC_USE_SHARED_LIBRARIES)
188: suffix: fgmres_geneo_20_p_2_geneo
189: output_file: output/ex76_fgmres_geneo_20_p_2.out
190: nsize: 4
191: 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} -mat_type {{aij sbaij}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
192: # PCHPDDM + KSPHPDDM test to exercise multilevel + multiple RHS in one go
193: test:
194: requires: hpddm slepc datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) define(PETSC_HAVE_DYNAMIC_LIBRARIES) define(PETSC_USE_SHARED_LIBRARIES)
195: suffix: fgmres_geneo_20_p_2_geneo_rhs
196: output_file: output/ex76_fgmres_geneo_20_p_2.out
197: # for -pc_hpddm_coarse_correction additive
198: filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 37/Linear solve converged due to CONVERGED_RTOL iterations 25/g"
199: nsize: 4
200: 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}
202: TEST*/