1 static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
2 Input arguments are:\n\
3 -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";
4
5 #include <petscksp.h>
6 #include <petsclog.h>
7
KSPTestResidualMonitor(KSP ksp,PetscInt i,PetscReal r,PetscCtx ctx)8 static PetscErrorCode KSPTestResidualMonitor(KSP ksp, PetscInt i, PetscReal r, PetscCtx ctx)
9 {
10 Vec *t, *v;
11 PetscReal err;
12
13 PetscFunctionBeginUser;
14 PetscCall(KSPCreateVecs(ksp, 2, &t, 2, &v));
15 PetscCall(KSPBuildResidualDefault(ksp, t[0], v[0], &v[0]));
16 PetscCall(KSPBuildResidual(ksp, t[1], v[1], &v[1]));
17 PetscCall(VecAXPY(v[1], -1.0, v[0]));
18 PetscCall(VecNorm(v[1], NORM_INFINITY, &err));
19 PetscCheck(err <= PETSC_SMALL, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Inconsistent residual computed at step %" PetscInt_FMT ": %g (KSP %g)", i, (double)err, (double)r);
20 PetscCall(VecDestroyVecs(2, &t));
21 PetscCall(VecDestroyVecs(2, &v));
22 PetscFunctionReturn(PETSC_SUCCESS);
23 }
24
main(int argc,char ** args)25 int main(int argc, char **args)
26 {
27 PetscInt its;
28 PetscLogStage stage1, stage2;
29 PetscReal norm;
30 Vec x, b, u;
31 Mat A;
32 char file[PETSC_MAX_PATH_LEN];
33 PetscViewer fd;
34 PetscBool table = PETSC_FALSE, flg, test_residual = PETSC_FALSE, b_in_f = PETSC_TRUE;
35 KSP ksp;
36
37 PetscFunctionBeginUser;
38 PetscCall(PetscInitialize(&argc, &args, NULL, help));
39 PetscCall(PetscOptionsGetBool(NULL, NULL, "-table", &table, NULL));
40 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_residual", &test_residual, NULL));
41 PetscCall(PetscOptionsGetBool(NULL, NULL, "-b_in_f", &b_in_f, NULL));
42
43 /* Read matrix and RHS */
44 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
45 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must indicate binary file with the -f option");
46 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
47 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
48 PetscCall(MatLoad(A, fd));
49 if (b_in_f) {
50 PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
51 PetscCall(VecLoad(b, fd));
52 } else {
53 PetscCall(MatCreateVecs(A, NULL, &b));
54 PetscCall(VecSetRandom(b, NULL));
55 }
56 PetscCall(PetscViewerDestroy(&fd));
57
58 /*
59 If the load matrix is larger than the vector, due to being padded
60 to match the blocksize then create a new padded vector
61 */
62 {
63 PetscInt m, n, j, mvec, start, end, indx;
64 Vec tmp;
65 PetscScalar *bold;
66
67 PetscCall(MatGetLocalSize(A, &m, &n));
68 PetscCall(VecCreate(PETSC_COMM_WORLD, &tmp));
69 PetscCall(VecSetSizes(tmp, m, PETSC_DECIDE));
70 PetscCall(VecSetFromOptions(tmp));
71 PetscCall(VecGetOwnershipRange(b, &start, &end));
72 PetscCall(VecGetLocalSize(b, &mvec));
73 PetscCall(VecGetArray(b, &bold));
74 for (j = 0; j < mvec; j++) {
75 indx = start + j;
76 PetscCall(VecSetValues(tmp, 1, &indx, bold + j, INSERT_VALUES));
77 }
78 PetscCall(VecRestoreArray(b, &bold));
79 PetscCall(VecDestroy(&b));
80 PetscCall(VecAssemblyBegin(tmp));
81 PetscCall(VecAssemblyEnd(tmp));
82 b = tmp;
83 }
84 PetscCall(VecDuplicate(b, &x));
85 PetscCall(VecDuplicate(b, &u));
86
87 PetscCall(VecSet(x, 0.0));
88 PetscCall(PetscBarrier((PetscObject)A));
89
90 PetscCall(PetscLogStageRegister("mystage 1", &stage1));
91 PetscCall(PetscLogStagePush(stage1));
92 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
93 PetscCall(KSPSetOperators(ksp, A, A));
94 PetscCall(KSPSetFromOptions(ksp));
95 if (test_residual) PetscCall(KSPMonitorSet(ksp, KSPTestResidualMonitor, NULL, NULL));
96 PetscCall(KSPSetUp(ksp));
97 PetscCall(KSPSetUpOnBlocks(ksp));
98 PetscCall(PetscLogStagePop());
99 PetscCall(PetscBarrier((PetscObject)A));
100
101 PetscCall(PetscLogStageRegister("mystage 2", &stage2));
102 PetscCall(PetscLogStagePush(stage2));
103 PetscCall(KSPSolve(ksp, b, x));
104 PetscCall(PetscLogStagePop());
105
106 /* Show result */
107 PetscCall(MatMult(A, x, u));
108 PetscCall(VecAXPY(u, -1.0, b));
109 PetscCall(VecNorm(u, NORM_2, &norm));
110 PetscCall(KSPGetIterationNumber(ksp, &its));
111 /* matrix PC KSP Options its residual */
112 if (table) {
113 char *matrixname, kspinfo[120];
114 PetscViewer viewer;
115 PetscCall(PetscViewerStringOpen(PETSC_COMM_WORLD, kspinfo, sizeof(kspinfo), &viewer));
116 PetscCall(KSPView(ksp, viewer));
117 PetscCall(PetscStrrchr(file, '/', &matrixname));
118 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-8.8s %3" PetscInt_FMT " %2.0e %s \n", matrixname, its, (double)norm, kspinfo));
119 PetscCall(PetscViewerDestroy(&viewer));
120 } else {
121 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
122 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm = %g\n", (double)norm));
123 }
124
125 /* Cleanup */
126 PetscCall(KSPDestroy(&ksp));
127 PetscCall(VecDestroy(&x));
128 PetscCall(VecDestroy(&b));
129 PetscCall(VecDestroy(&u));
130 PetscCall(MatDestroy(&A));
131 PetscCall(PetscFinalize());
132 return 0;
133 }
134
135 /*TEST
136
137 test:
138 args: -ksp_type preonly -pc_type lu -options_left no -f ${DATAFILESPATH}/matrices/arco1
139 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
140
141 test:
142 suffix: 2
143 args: -sub_pc_type ilu -options_left no -f ${DATAFILESPATH}/matrices/arco1 -ksp_gmres_restart 100 -ksp_gmres_cgs_refinement_type refine_always -sub_ksp_type preonly -pc_type bjacobi -pc_bjacobi_blocks 8 -sub_pc_factor_in_place -ksp_monitor_short
144 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
145
146 test:
147 suffix: 7
148 args: -ksp_gmres_cgs_refinement_type refine_always -pc_type asm -pc_asm_blocks 6 -f ${DATAFILESPATH}/matrices/small -matload_block_size 6 -ksp_monitor_short
149 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
150
151 test:
152 requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
153 suffix: 3
154 filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
155 args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -pc_type none -ksp_type {{cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres pipefgmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs qmrcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr lsqr qcg bicg minres symmlq lcd gcr pipegcr cgls}} -ksp_max_it 20 -ksp_error_if_not_converged -ksp_converged_reason -test_residual
156
157 test:
158 requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
159 suffix: 3_maxits
160 output_file: output/ex6_maxits.out
161 args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -pc_type none -ksp_type {{chebyshev cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres pipefgmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs qmrcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr qcg bicg minres symmlq lcd gcr pipegcr cgls richardson}} -ksp_max_it 4 -ksp_error_if_not_converged -ksp_converged_maxits -ksp_converged_reason -test_residual -ksp_norm_type none
162
163 testset:
164 requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
165 output_file: output/ex6_skip.out
166 args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -pc_type none -ksp_max_it 8 -ksp_error_if_not_converged -ksp_convergence_test skip -ksp_converged_reason -test_residual
167 #SYMMLQ converges in 4 iterations and then generate nans
168 test:
169 suffix: 3_skip
170 args: -ksp_type {{chebyshev cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs qmrcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr qcg bicg minres lcd gcr cgls richardson}}
171 #PIPEGCR generates nans on linux-knl
172 #PIPEFGMRES can have happy breakdown which is not handled well with no convergence test
173 test:
174 requires: !defined(PETSC_USE_AVX512_KERNELS)
175 suffix: 3_skip_pipegcr
176 args: -ksp_type pipegcr
177 test:
178 requires: hpddm
179 suffix: 3_skip_hpddm
180 args: -ksp_type hpddm -ksp_hpddm_type {{cg gmres bgmres bcg bfbcg gcrodr bgcrodr}}
181
182 test:
183 requires: double !complex !defined(PETSC_USE_64BIT_INDICES) hpddm
184 suffix: 3_hpddm
185 output_file: output/ex6_3.out
186 filter: sed -e "s/CONVERGED_RTOL/CONVERGED_ATOL/g"
187 args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/spd-real-int32-float64 -pc_type none -ksp_type hpddm -ksp_hpddm_type {{cg gmres bgmres bcg bfbcg gcrodr bgcrodr}} -ksp_max_it 20 -ksp_error_if_not_converged -ksp_converged_reason -test_residual
188
189 # test CG shortcut for residual access
190 test:
191 suffix: 4
192 args: -ksp_converged_reason -ksp_max_it 20 -ksp_converged_maxits -ksp_type {{cg pipecg groppcg}} -ksp_norm_type {{preconditioned unpreconditioned natural}separate output} -pc_type {{bjacobi none}separate output} -f ${DATAFILESPATH}/matrices/poisson_2d13p -b_in_f 0 -test_residual
193 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
194
195 TEST*/
196