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: PetscViewer viewer;
18: char name[PETSC_MAX_PATH_LEN], type[256];
19: PetscBool breakdown = PETSC_FALSE, flg;
20: KSPConvergedReason reason;
22: PetscFunctionBeginUser;
23: PetscCall(PetscInitialize(&argc, &args, NULL, help));
24: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
25: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
26: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", name, sizeof(name), &flg));
27: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must provide a binary file for the matrix");
28: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
29: PetscCall(PetscOptionsGetBool(NULL, NULL, "-breakdown", &breakdown, NULL));
30: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ksp_hpddm_deflation_tol", &deflation, NULL));
31: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
32: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
33: PetscCall(KSPSetOperators(ksp, A, A));
34: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
35: PetscCall(MatLoad(A, viewer));
36: PetscCall(PetscViewerDestroy(&viewer));
37: PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
38: PetscCall(PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, MATAIJ, type, sizeof(type), &flg));
39: PetscOptionsEnd();
40: if (flg) {
41: PetscCall(PetscStrcmp(type, MATKAIJ, &flg));
42: if (!flg) {
43: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
44: PetscCall(MatConvert(A, type, MAT_INPLACE_MATRIX, &A));
45: } else {
46: if (size > 2) {
47: PetscCall(MatGetSize(A, &M, NULL));
48: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
49: if (rank > 1) PetscCall(MatSetSizes(B, 0, 0, M, M));
50: else PetscCall(MatSetSizes(B, rank ? M - M / 2 : M / 2, rank ? M - M / 2 : M / 2, M, M));
51: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
52: PetscCall(MatLoad(B, viewer));
53: PetscCall(PetscViewerDestroy(&viewer));
54: PetscCall(MatHeaderReplace(A, &B));
55: }
56: PetscCall(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: PetscCall(MatCreateKAIJ(A, N, N, S, T, &KA));
62: }
63: }
64: if (!flg) {
65: if (size > 4) {
66: Mat B;
67: PetscCall(MatGetSize(A, &M, NULL));
68: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
69: if (rank > 3) PetscCall(MatSetSizes(B, 0, 0, M, M));
70: else PetscCall(MatSetSizes(B, rank == 0 ? M - 3 * (M / 4) : M / 4, rank == 0 ? M - 3 * (M / 4) : M / 4, M, M));
71: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
72: PetscCall(MatLoad(B, viewer));
73: PetscCall(PetscViewerDestroy(&viewer));
74: PetscCall(MatHeaderReplace(A, &B));
75: }
76: }
77: PetscCall(MatGetLocalSize(A, &m, NULL));
78: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
79: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &X));
80: if (!breakdown) PetscCall(MatSetRandom(B, NULL));
81: PetscCall(KSPSetFromOptions(ksp));
82: if (!flg) {
83: if (!breakdown) {
84: PetscCall(KSPMatSolve(ksp, B, X));
85: PetscCall(KSPGetMatSolveBatchSize(ksp, &M));
86: if (M != PETSC_DECIDE) {
87: PetscCall(KSPSetMatSolveBatchSize(ksp, PETSC_DECIDE));
88: PetscCall(MatZeroEntries(X));
89: PetscCall(KSPMatSolve(ksp, B, X));
90: }
91: PetscCall(KSPGetPC(ksp, &pc));
92: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLU, &flg));
93: if (flg) {
94: PetscCall(PCFactorGetMatrix(pc, &F));
95: PetscCall(MatMatSolve(F, B, B));
96: PetscCall(MatAYPX(B, -1.0, X, SAME_NONZERO_PATTERN));
97: PetscCall(MatNorm(B, NORM_INFINITY, &norm));
98: PetscCheck(norm < 100 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPMatSolve() and MatMatSolve() difference has nonzero norm %g", (double)norm);
99: }
100: } else {
101: PetscCall(MatZeroEntries(B));
102: PetscCall(KSPMatSolve(ksp, B, X));
103: PetscCall(KSPGetConvergedReason(ksp, &reason));
104: PetscCheck(reason == KSP_CONVERGED_HAPPY_BREAKDOWN, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPConvergedReason() %s != KSP_CONVERGED_HAPPY_BREAKDOWN", KSPConvergedReasons[reason]);
105: PetscCall(MatDenseGetArrayWrite(B, &x));
106: for (i = 0; i < m * N; ++i) x[i] = 1.0;
107: PetscCall(MatDenseRestoreArrayWrite(B, &x));
108: PetscCall(KSPMatSolve(ksp, B, X));
109: PetscCall(KSPGetConvergedReason(ksp, &reason));
110: PetscCheck(reason == KSP_DIVERGED_BREAKDOWN || deflation >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPConvergedReason() %s != KSP_DIVERGED_BREAKDOWN", KSPConvergedReasons[reason]);
111: }
112: } else {
113: PetscCall(KSPSetOperators(ksp, KA, KA));
114: PetscCall(MatGetSize(KA, &M, NULL));
115: PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, m * N, M, &cb));
116: PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, m * N, M, &cx));
117: PetscCall(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: PetscCall(KSPSolve(ksp, cb, cx));
121: PetscCall(VecDestroy(&cx));
122: PetscCall(VecDestroy(&cb));
123: }
124: PetscCall(MatDestroy(&X));
125: PetscCall(MatDestroy(&B));
126: PetscCall(PetscFree2(S, T));
127: PetscCall(MatDestroy(&KA));
128: PetscCall(MatDestroy(&A));
129: PetscCall(KSPDestroy(&ksp));
130: PetscCall(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
162: output_file: output/empty.out
164: testset:
165: nsize: 1
166: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
167: args: -N 3 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_type hpddm -breakdown
168: test:
169: suffix: breakdown_wo_deflation
170: output_file: output/empty.out
171: args: -pc_type none -ksp_hpddm_type {{bcg bgmres bgcrodr bfbcg}shared output}
172: test:
173: suffix: breakdown_w_deflation
174: output_file: output/ex77_deflation.out
175: filter: sed "s/BGCRODR/BGMRES/g"
176: args: -pc_type lu -ksp_hpddm_type {{bgmres bgcrodr}shared output} -ksp_hpddm_deflation_tol 1e-1 -info :ksp
178: test:
179: nsize: 2
180: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
181: 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}
183: test:
184: nsize: 3
185: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
186: suffix: kaij_zero
187: output_file: output/ex77_ksp_hpddm_type-bgmres.out
188: 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
190: test:
191: nsize: 4
192: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
193: suffix: 4_slepc
194: output_file: output/ex77_4.out
195: filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
196: 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
198: testset:
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: filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
202: 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
203: test:
204: requires: elemental
205: suffix: 4_elemental
206: output_file: output/ex77_4.out
207: args: -ksp_hpddm_recycle_mat_type elemental
208: test:
209: requires: elemental
210: suffix: 4_elemental_matmat
211: output_file: output/ex77_4.out
212: 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
213: test:
214: requires: scalapack
215: suffix: 4_scalapack
216: output_file: output/ex77_4.out
217: args: -ksp_hpddm_recycle_mat_type scalapack
219: test:
220: nsize: 5
221: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
222: suffix: 4_zero
223: output_file: output/ex77_4.out
224: 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
226: TEST*/