Actual source code: ex77.c
1: #include <petsc.h>
3: static char help[] = "Solves a linear system with a block of right-hand sides using KSPHPDDM.\n\n";
5: int main(int argc,char **args)
6: {
7: Mat X,B; /* computed solutions and RHS */
8: Vec cx,cb; /* columns of X and B */
9: Mat A,KA = NULL; /* linear system matrix */
10: KSP ksp; /* linear solver context */
11: PC pc; /* preconditioner context */
12: Mat F; /* factored matrix from the preconditioner context */
13: PetscScalar *x,*S = NULL,*T = NULL;
14: PetscReal norm,deflation = -1.0;
15: PetscInt m,M,N = 5,i;
16: PetscMPIInt rank,size;
17: const char *deft = MATAIJ;
18: PetscViewer viewer;
19: char name[PETSC_MAX_PATH_LEN],type[256];
20: PetscBool breakdown = PETSC_FALSE,flg;
21: KSPConvergedReason reason;
23: PetscInitialize(&argc,&args,NULL,help);
24: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
25: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: PetscOptionsGetString(NULL,NULL,"-f",name,sizeof(name),&flg);
28: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
29: PetscOptionsGetBool(NULL,NULL,"-breakdown",&breakdown,NULL);
30: PetscOptionsGetReal(NULL,NULL,"-ksp_hpddm_deflation_tol",&deflation,NULL);
31: MatCreate(PETSC_COMM_WORLD,&A);
32: KSPCreate(PETSC_COMM_WORLD,&ksp);
33: KSPSetOperators(ksp,A,A);
34: PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);
35: MatLoad(A,viewer);
36: PetscViewerDestroy(&viewer);
37: PetscOptionsBegin(PETSC_COMM_WORLD,"","","");
38: PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
39: PetscOptionsEnd();
40: if (flg) {
41: PetscStrcmp(type,MATKAIJ,&flg);
42: if (!flg) {
43: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
44: MatConvert(A,type,MAT_INPLACE_MATRIX,&A);
45: } else {
46: if (size > 2) {
47: MatGetSize(A,&M,NULL);
48: MatCreate(PETSC_COMM_WORLD,&B);
49: if (rank > 1) MatSetSizes(B,0,0,M,M);
50: else MatSetSizes(B,rank?M-M/2:M/2,rank?M-M/2:M/2,M,M);
51: PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);
52: MatLoad(B,viewer);
53: PetscViewerDestroy(&viewer);
54: MatHeaderReplace(A,&B);
55: }
56: PetscCalloc2(N*N,&S,N*N,&T);
57: for (i=0; i<N; i++) { /* really easy problem used for testing */
58: S[i*(N+1)] = 1e+6;
59: T[i*(N+1)] = 1e-2;
60: }
61: MatCreateKAIJ(A,N,N,S,T,&KA);
62: }
63: }
64: if (!flg) {
65: if (size > 4) {
66: Mat B;
67: MatGetSize(A,&M,NULL);
68: MatCreate(PETSC_COMM_WORLD,&B);
69: if (rank > 3) MatSetSizes(B,0,0,M,M);
70: else MatSetSizes(B,rank == 0?M-3*(M/4):M/4,rank == 0?M-3*(M/4):M/4,M,M);
71: PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);
72: MatLoad(B,viewer);
73: PetscViewerDestroy(&viewer);
74: MatHeaderReplace(A,&B);
75: }
76: }
77: MatGetLocalSize(A,&m,NULL);
78: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&B);
79: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&X);
80: if (!breakdown) MatSetRandom(B,NULL);
81: KSPSetFromOptions(ksp);
82: if (!flg) {
83: if (!breakdown) {
84: KSPMatSolve(ksp,B,X);
85: KSPGetMatSolveBatchSize(ksp,&M);
86: if (M != PETSC_DECIDE) {
87: KSPSetMatSolveBatchSize(ksp,PETSC_DECIDE);
88: MatZeroEntries(X);
89: KSPMatSolve(ksp,B,X);
90: }
91: KSPGetPC(ksp,&pc);
92: PetscObjectTypeCompare((PetscObject)pc,PCLU,&flg);
93: if (flg) {
94: PCFactorGetMatrix(pc,&F);
95: MatMatSolve(F,B,B);
96: MatAYPX(B,-1.0,X,SAME_NONZERO_PATTERN);
97: MatNorm(B,NORM_INFINITY,&norm);
99: }
100: } else {
101: MatZeroEntries(B);
102: KSPMatSolve(ksp,B,X);
103: KSPGetConvergedReason(ksp,&reason);
105: MatDenseGetArrayWrite(B,&x);
106: for (i=0; i<m*N; ++i) x[i] = 1.0;
107: MatDenseRestoreArrayWrite(B,&x);
108: KSPMatSolve(ksp,B,X);
109: KSPGetConvergedReason(ksp,&reason);
111: }
112: } else {
113: KSPSetOperators(ksp,KA,KA);
114: MatGetSize(KA,&M,NULL);
115: VecCreateMPI(PETSC_COMM_WORLD,m*N,M,&cb);
116: VecCreateMPI(PETSC_COMM_WORLD,m*N,M,&cx);
117: VecSetRandom(cb,NULL);
118: /* solving with MatKAIJ is equivalent to block solving with row-major RHS and solutions */
119: /* only applies if MatKAIJGetScaledIdentity() returns true */
120: KSPSolve(ksp,cb,cx);
121: VecDestroy(&cx);
122: VecDestroy(&cb);
123: }
124: MatDestroy(&X);
125: MatDestroy(&B);
126: PetscFree2(S,T);
127: MatDestroy(&KA);
128: MatDestroy(&A);
129: KSPDestroy(&ksp);
130: PetscFinalize();
131: return 0;
132: }
134: /*TEST
136: testset:
137: nsize: 2
138: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
139: args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -mat_type {{aij sbaij}shared output}
140: test:
141: suffix: 1
142: args:
143: test:
144: suffix: 2
145: requires: hpddm
146: args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type {{gmres bgmres}separate output}
147: test:
148: suffix: 3
149: requires: hpddm
150: args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type {{gcrodr bgcrodr}separate output}
151: test:
152: nsize: 4
153: suffix: 4
154: requires: hpddm
155: args: -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_batch_size 5
157: test:
158: nsize: 1
159: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
160: suffix: preonly
161: args: -N 6 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -pc_type lu -ksp_type hpddm -ksp_hpddm_type preonly
163: testset:
164: nsize: 1
165: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
166: args: -N 3 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_type hpddm -breakdown
167: test:
168: suffix: breakdown_wo_deflation
169: output_file: output/ex77_preonly.out
170: args: -pc_type none -ksp_hpddm_type {{bcg bgmres bgcrodr bfbcg}shared output}
171: test:
172: suffix: breakdown_w_deflation
173: output_file: output/ex77_deflation.out
174: filter: sed "s/BGCRODR/BGMRES/g"
175: args: -pc_type lu -ksp_hpddm_type {{bgmres bgcrodr}shared output} -ksp_hpddm_deflation_tol 1e-1 -info :ksp
177: test:
178: nsize: 2
179: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
180: args: -N 12 -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -mat_type kaij -pc_type pbjacobi -ksp_type hpddm -ksp_hpddm_type {{gmres bgmres}separate output}
182: test:
183: nsize: 3
184: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
185: suffix: kaij_zero
186: output_file: output/ex77_ksp_hpddm_type-bgmres.out
187: args: -N 12 -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -mat_type kaij -pc_type pbjacobi -ksp_type hpddm -ksp_hpddm_type bgmres
189: test:
190: nsize: 4
191: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
192: suffix: 4_slepc
193: output_file: output/ex77_4.out
194: filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
195: args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_batch_size 5 -ksp_hpddm_recycle_redistribute 2 -ksp_hpddm_recycle_mat_type {{aij dense}shared output} -ksp_hpddm_recycle_eps_converged_reason -ksp_hpddm_recycle_st_pc_type redundant
197: testset:
198: nsize: 4
199: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
200: filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
201: args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_batch_size 5 -ksp_hpddm_recycle_redistribute 2 -ksp_hpddm_recycle_eps_converged_reason
202: test:
203: requires: elemental
204: suffix: 4_elemental
205: output_file: output/ex77_4.out
206: args: -ksp_hpddm_recycle_mat_type elemental
207: test:
208: requires: elemental
209: suffix: 4_elemental_matmat
210: output_file: output/ex77_4.out
211: args: -ksp_hpddm_recycle_mat_type elemental -ksp_hpddm_recycle_eps_type subspace -ksp_hpddm_recycle_eps_target 0 -ksp_hpddm_recycle_eps_target_magnitude -ksp_hpddm_recycle_st_type sinvert -ksp_hpddm_recycle_bv_type mat -ksp_hpddm_recycle_bv_orthog_block svqb
212: test:
213: requires: scalapack
214: suffix: 4_scalapack
215: output_file: output/ex77_4.out
216: args: -ksp_hpddm_recycle_mat_type scalapack
218: test:
219: nsize: 5
220: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
221: suffix: 4_zero
222: output_file: output/ex77_4.out
223: args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_batch_size 5
225: TEST*/