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: PetscFunctionBeginUser;
24: PetscCall(PetscInitialize(&argc, &args, NULL, help));
25: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
26: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
27: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", name, sizeof(name), &flg));
28: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must provide a binary file for the matrix");
29: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
30: PetscCall(PetscOptionsGetBool(NULL, NULL, "-breakdown", &breakdown, NULL));
31: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ksp_hpddm_deflation_tol", &deflation, NULL));
32: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
33: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
34: PetscCall(KSPSetOperators(ksp, A, A));
35: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
36: PetscCall(MatLoad(A, viewer));
37: PetscCall(PetscViewerDestroy(&viewer));
38: PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
39: PetscCall(PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, deft, type, 256, &flg));
40: PetscOptionsEnd();
41: if (flg) {
42: PetscCall(PetscStrcmp(type, MATKAIJ, &flg));
43: if (!flg) {
44: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
45: PetscCall(MatConvert(A, type, MAT_INPLACE_MATRIX, &A));
46: } else {
47: if (size > 2) {
48: PetscCall(MatGetSize(A, &M, NULL));
49: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
50: if (rank > 1) PetscCall(MatSetSizes(B, 0, 0, M, M));
51: else PetscCall(MatSetSizes(B, rank ? M - M / 2 : M / 2, rank ? M - M / 2 : M / 2, M, M));
52: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
53: PetscCall(MatLoad(B, viewer));
54: PetscCall(PetscViewerDestroy(&viewer));
55: PetscCall(MatHeaderReplace(A, &B));
56: }
57: PetscCall(PetscCalloc2(N * N, &S, N * N, &T));
58: for (i = 0; i < N; i++) { /* really easy problem used for testing */
59: S[i * (N + 1)] = 1e+6;
60: T[i * (N + 1)] = 1e-2;
61: }
62: PetscCall(MatCreateKAIJ(A, N, N, S, T, &KA));
63: }
64: }
65: if (!flg) {
66: if (size > 4) {
67: Mat B;
68: PetscCall(MatGetSize(A, &M, NULL));
69: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
70: if (rank > 3) PetscCall(MatSetSizes(B, 0, 0, M, M));
71: else PetscCall(MatSetSizes(B, rank == 0 ? M - 3 * (M / 4) : M / 4, rank == 0 ? M - 3 * (M / 4) : M / 4, M, M));
72: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
73: PetscCall(MatLoad(B, viewer));
74: PetscCall(PetscViewerDestroy(&viewer));
75: PetscCall(MatHeaderReplace(A, &B));
76: }
77: }
78: PetscCall(MatGetLocalSize(A, &m, NULL));
79: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
80: PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &X));
81: if (!breakdown) PetscCall(MatSetRandom(B, NULL));
82: PetscCall(KSPSetFromOptions(ksp));
83: if (!flg) {
84: if (!breakdown) {
85: PetscCall(KSPMatSolve(ksp, B, X));
86: PetscCall(KSPGetMatSolveBatchSize(ksp, &M));
87: if (M != PETSC_DECIDE) {
88: PetscCall(KSPSetMatSolveBatchSize(ksp, PETSC_DECIDE));
89: PetscCall(MatZeroEntries(X));
90: PetscCall(KSPMatSolve(ksp, B, X));
91: }
92: PetscCall(KSPGetPC(ksp, &pc));
93: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLU, &flg));
94: if (flg) {
95: PetscCall(PCFactorGetMatrix(pc, &F));
96: PetscCall(MatMatSolve(F, B, B));
97: PetscCall(MatAYPX(B, -1.0, X, SAME_NONZERO_PATTERN));
98: PetscCall(MatNorm(B, NORM_INFINITY, &norm));
99: PetscCheck(norm < 100 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPMatSolve() and MatMatSolve() difference has nonzero norm %g", (double)norm);
100: }
101: } else {
102: PetscCall(MatZeroEntries(B));
103: PetscCall(KSPMatSolve(ksp, B, X));
104: PetscCall(KSPGetConvergedReason(ksp, &reason));
105: PetscCheck(reason == KSP_CONVERGED_HAPPY_BREAKDOWN, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPConvergedReason() %s != KSP_CONVERGED_HAPPY_BREAKDOWN", KSPConvergedReasons[reason]);
106: PetscCall(MatDenseGetArrayWrite(B, &x));
107: for (i = 0; i < m * N; ++i) x[i] = 1.0;
108: PetscCall(MatDenseRestoreArrayWrite(B, &x));
109: PetscCall(KSPMatSolve(ksp, B, X));
110: PetscCall(KSPGetConvergedReason(ksp, &reason));
111: PetscCheck(reason == KSP_DIVERGED_BREAKDOWN || deflation >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPConvergedReason() %s != KSP_DIVERGED_BREAKDOWN", KSPConvergedReasons[reason]);
112: }
113: } else {
114: PetscCall(KSPSetOperators(ksp, KA, KA));
115: PetscCall(MatGetSize(KA, &M, NULL));
116: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, m * N, M, &cb));
117: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, m * N, M, &cx));
118: PetscCall(VecSetRandom(cb, NULL));
119: /* solving with MatKAIJ is equivalent to block solving with row-major RHS and solutions */
120: /* only applies if MatKAIJGetScaledIdentity() returns true */
121: PetscCall(KSPSolve(ksp, cb, cx));
122: PetscCall(VecDestroy(&cx));
123: PetscCall(VecDestroy(&cb));
124: }
125: PetscCall(MatDestroy(&X));
126: PetscCall(MatDestroy(&B));
127: PetscCall(PetscFree2(S, T));
128: PetscCall(MatDestroy(&KA));
129: PetscCall(MatDestroy(&A));
130: PetscCall(KSPDestroy(&ksp));
131: PetscCall(PetscFinalize());
132: return 0;
133: }
135: /*TEST
137: testset:
138: nsize: 2
139: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
140: args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -mat_type {{aij sbaij}shared output}
141: test:
142: suffix: 1
143: args:
144: test:
145: suffix: 2
146: requires: hpddm
147: args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type {{gmres bgmres}separate output}
148: test:
149: suffix: 3
150: requires: hpddm
151: args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type {{gcrodr bgcrodr}separate output}
152: test:
153: nsize: 4
154: suffix: 4
155: requires: hpddm
156: 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
158: test:
159: nsize: 1
160: requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
161: suffix: preonly
162: args: -N 6 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -pc_type lu -ksp_type hpddm -ksp_hpddm_type preonly
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/ex77_preonly.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*/