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;
22: PetscErrorCode ierr;
24: PetscInitialize(&argc,&args,NULL,help);
25: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
26: MPI_Comm_size(PETSC_COMM_WORLD,&size);
27: PetscOptionsGetString(NULL,NULL,"-f",name,sizeof(name),&flg);
29: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
30: PetscOptionsGetBool(NULL,NULL,"-breakdown",&breakdown,NULL);
31: PetscOptionsGetReal(NULL,NULL,"-ksp_hpddm_deflation_tol",&deflation,NULL);
32: MatCreate(PETSC_COMM_WORLD,&A);
33: KSPCreate(PETSC_COMM_WORLD,&ksp);
34: KSPSetOperators(ksp,A,A);
35: PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);
36: MatLoad(A,viewer);
37: PetscViewerDestroy(&viewer);
38: PetscOptionsBegin(PETSC_COMM_WORLD,"","","");
39: PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,type,256,&flg);
40: PetscOptionsEnd();
41: if (flg) {
42: PetscStrcmp(type,MATKAIJ,&flg);
43: if (!flg) {
44: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
45: MatConvert(A,type,MAT_INPLACE_MATRIX,&A);
46: } else {
47: if (size > 2) {
48: MatGetSize(A,&M,NULL);
49: MatCreate(PETSC_COMM_WORLD,&B);
50: if (rank > 1) {
51: MatSetSizes(B,0,0,M,M);
52: } else {
53: MatSetSizes(B,rank?M-M/2:M/2,rank?M-M/2:M/2,M,M);
54: }
55: PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);
56: MatLoad(B,viewer);
57: PetscViewerDestroy(&viewer);
58: MatHeaderReplace(A,&B);
59: }
60: PetscCalloc2(N*N,&S,N*N,&T);
61: for (i=0; i<N; i++) { /* really easy problem used for testing */
62: S[i*(N+1)] = 1e+6;
63: T[i*(N+1)] = 1e-2;
64: }
65: MatCreateKAIJ(A,N,N,S,T,&KA);
66: }
67: }
68: if (!flg) {
69: if (size > 4) {
70: Mat B;
71: MatGetSize(A,&M,NULL);
72: MatCreate(PETSC_COMM_WORLD,&B);
73: if (rank > 3) {
74: MatSetSizes(B,0,0,M,M);
75: } else {
76: MatSetSizes(B,rank == 0?M-3*(M/4):M/4,rank == 0?M-3*(M/4):M/4,M,M);
77: }
78: PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,&viewer);
79: MatLoad(B,viewer);
80: PetscViewerDestroy(&viewer);
81: MatHeaderReplace(A,&B);
82: }
83: }
84: MatGetLocalSize(A,&m,NULL);
85: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&B);
86: MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,NULL,&X);
87: if (!breakdown) {
88: MatSetRandom(B,NULL);
89: }
90: KSPSetFromOptions(ksp);
91: if (!flg) {
92: if (!breakdown) {
93: KSPMatSolve(ksp,B,X);
94: KSPGetMatSolveBatchSize(ksp,&M);
95: if (M != PETSC_DECIDE) {
96: KSPSetMatSolveBatchSize(ksp,PETSC_DECIDE);
97: MatZeroEntries(X);
98: KSPMatSolve(ksp,B,X);
99: }
100: KSPGetPC(ksp,&pc);
101: PetscObjectTypeCompare((PetscObject)pc,PCLU,&flg);
102: if (flg) {
103: PCFactorGetMatrix(pc,&F);
104: MatMatSolve(F,B,B);
105: MatAYPX(B,-1.0,X,SAME_NONZERO_PATTERN);
106: MatNorm(B,NORM_INFINITY,&norm);
108: }
109: } else {
110: MatZeroEntries(B);
111: KSPMatSolve(ksp,B,X);
112: KSPGetConvergedReason(ksp,&reason);
114: MatDenseGetArrayWrite(B,&x);
115: for (i=0; i<m*N; ++i) x[i] = 1.0;
116: MatDenseRestoreArrayWrite(B,&x);
117: KSPMatSolve(ksp,B,X);
118: KSPGetConvergedReason(ksp,&reason);
120: }
121: } else {
122: KSPSetOperators(ksp,KA,KA);
123: MatGetSize(KA,&M,NULL);
124: VecCreateMPI(PETSC_COMM_WORLD,m*N,M,&cb);
125: VecCreateMPI(PETSC_COMM_WORLD,m*N,M,&cx);
126: VecSetRandom(cb,NULL);
127: /* solving with MatKAIJ is equivalent to block solving with row-major RHS and solutions */
128: /* only applies if MatKAIJGetScaledIdentity() returns true */
129: KSPSolve(ksp,cb,cx);
130: VecDestroy(&cx);
131: VecDestroy(&cb);
132: }
133: MatDestroy(&X);
134: MatDestroy(&B);
135: PetscFree2(S,T);
136: MatDestroy(&KA);
137: MatDestroy(&A);
138: KSPDestroy(&ksp);
139: PetscFinalize();
140: return 0;
141: }
143: /*TEST
145: testset:
146: nsize: 2
147: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
148: args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -mat_type {{aij sbaij}shared output}
149: test:
150: suffix: 1
151: args:
152: test:
153: suffix: 2
154: requires: hpddm
155: args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type {{gmres bgmres}separate output}
156: test:
157: suffix: 3
158: requires: hpddm
159: args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type {{gcrodr bgcrodr}separate output}
160: test:
161: nsize: 4
162: suffix: 4
163: requires: hpddm
164: 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
166: test:
167: nsize: 1
168: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
169: suffix: preonly
170: args: -N 6 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -pc_type lu -ksp_type hpddm -ksp_hpddm_type preonly
172: testset:
173: nsize: 1
174: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
175: args: -N 3 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_type hpddm -breakdown
176: test:
177: suffix: breakdown_wo_deflation
178: output_file: output/ex77_preonly.out
179: args: -pc_type none -ksp_hpddm_type {{bcg bgmres bgcrodr bfbcg}shared output}
180: test:
181: suffix: breakdown_w_deflation
182: output_file: output/ex77_deflation.out
183: filter: sed "s/BGCRODR/BGMRES/g"
184: args: -pc_type lu -ksp_hpddm_type {{bgmres bgcrodr}shared output} -ksp_hpddm_deflation_tol 1e-1 -info :ksp
186: test:
187: nsize: 2
188: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
189: 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}
191: test:
192: nsize: 3
193: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
194: suffix: kaij_zero
195: output_file: output/ex77_ksp_hpddm_type-bgmres.out
196: 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
198: test:
199: nsize: 4
200: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
201: suffix: 4_slepc
202: output_file: output/ex77_4.out
203: filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
204: 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
206: testset:
207: nsize: 4
208: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
209: filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
210: 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
211: test:
212: requires: elemental
213: suffix: 4_elemental
214: output_file: output/ex77_4.out
215: args: -ksp_hpddm_recycle_mat_type elemental
216: test:
217: requires: elemental
218: suffix: 4_elemental_matmat
219: output_file: output/ex77_4.out
220: 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
221: test:
222: requires: scalapack
223: suffix: 4_scalapack
224: output_file: output/ex77_4.out
225: args: -ksp_hpddm_recycle_mat_type scalapack
227: test:
228: nsize: 5
229: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
230: suffix: 4_zero
231: output_file: output/ex77_4.out
232: 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
234: TEST*/