1 static const char help[] = "Example demonstrating the effect of choosing FCG over CG\n\
2 for a simple diagonal system with a noisy preconditioner implemented using PCShell\n\
3 Accepts an option -n for the problem size\n\
4 Accepts an option -eta for the noise amplitude (set to 0 to deactivate)\n\
5 Accepts an option -diagfunc [1,2,3] to select from different eigenvalue distributions\n\
6 \n";
7
8 /*
9 Solve (in parallel) a diagonal linear system.
10
11 Using PCShell, we define a preconditioner which simply adds noise to the residual.
12
13 This example can be used to test the robustness of Krylov methods (particularly "flexible" ones for unitarily diagonalizable systems)
14 to varying preconditioners. Use the command line option -diagfunc [1,2,3] to choose between some predefined eigenvalue distributions.
15
16 The default behavior is to use a composite PC which combines (additively) an identity preconditioner with a preconditioner which
17 replaces the input with scaled noise.
18
19 To test with an inner Krylov method instead of noise, use PCKSP, e.g.
20 mpiexec -n 2 ./ex60 -eta 0 -ksp_type fcg -pc_type ksp -ksp_ksp_rtol 1e-1 -ksp_ksp_type cg -ksp_pc_type none
21 (note that eta is ignored here, and we specify the analogous quantity, the tolerance of the inner KSP solve,with -ksp_ksp_rtol)
22
23 To test by adding noise to a PC of your choosing (say ilu), run e.g.
24 mpiexec -n 2 ./ex60 -eta 0.1 -ksp_type fcg -sub_0_pc_type ilu
25
26 Contributed by Patrick Sanan
27 */
28
29 #include <petscksp.h>
30
31 /* Context to use with our noise PC */
32 typedef struct {
33 PetscReal eta;
34 PetscRandom random;
35 } PCNoise_Ctx;
36
PCApply_Noise(PC pc,Vec xin,Vec xout)37 PetscErrorCode PCApply_Noise(PC pc, Vec xin, Vec xout)
38 {
39 PCNoise_Ctx *ctx;
40 PetscReal nrmin, nrmnoise;
41
42 PetscFunctionBeginUser;
43 PetscCall(PCShellGetContext(pc, &ctx));
44
45 /* xout is ||xin|| * ctx->eta* f, where f is a pseudorandom unit vector
46 (Note that this should always be combined additively with another PC) */
47 PetscCall(VecSetRandom(xout, ctx->random));
48 PetscCall(VecNorm(xin, NORM_2, &nrmin));
49 PetscCall(VecNorm(xout, NORM_2, &nrmnoise));
50 PetscCall(VecScale(xout, ctx->eta * (nrmin / nrmnoise)));
51 PetscFunctionReturn(PETSC_SUCCESS);
52 }
53
PCSetup_Noise(PC pc)54 PetscErrorCode PCSetup_Noise(PC pc)
55 {
56 PCNoise_Ctx *ctx;
57
58 PetscFunctionBeginUser;
59 PetscCall(PCShellGetContext(pc, &ctx));
60 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &ctx->random));
61 PetscCall(PetscRandomSetInterval(ctx->random, -1.0, 1.0));
62 PetscFunctionReturn(PETSC_SUCCESS);
63 }
64
PCDestroy_Noise(PC pc)65 PetscErrorCode PCDestroy_Noise(PC pc)
66 {
67 PCNoise_Ctx *ctx;
68
69 PetscFunctionBeginUser;
70 PetscCall(PCShellGetContext(pc, &ctx));
71 PetscCall(PetscRandomDestroy(&ctx->random));
72 PetscFunctionReturn(PETSC_SUCCESS);
73 }
74
diagFunc1(PetscInt i,PetscInt n)75 PetscScalar diagFunc1(PetscInt i, PetscInt n)
76 {
77 const PetscScalar kappa = 5.0;
78 return 1.0 + (kappa * (PetscScalar)i) / (PetscScalar)(n - 1);
79 }
80
diagFunc2(PetscInt i,PetscInt n)81 PetscScalar diagFunc2(PetscInt i, PetscInt n)
82 {
83 const PetscScalar kappa = 50.0;
84 return 1.0 + (kappa * (PetscScalar)i) / (PetscScalar)(n - 1);
85 }
86
diagFunc3(PetscInt i,PetscInt n)87 PetscScalar diagFunc3(PetscInt i, PetscInt n)
88 {
89 const PetscScalar kappa = 10.0;
90 if (!i) {
91 return 1e-2;
92 } else {
93 return 1. + (kappa * ((PetscScalar)(i - 1))) / (PetscScalar)(n - 2);
94 }
95 }
96
AssembleDiagonalMatrix(Mat A,PetscScalar (* diagfunc)(PetscInt,PetscInt))97 static PetscErrorCode AssembleDiagonalMatrix(Mat A, PetscScalar (*diagfunc)(PetscInt, PetscInt))
98 {
99 PetscInt i, rstart, rend, n;
100 PetscScalar val;
101
102 PetscFunctionBeginUser;
103 PetscCall(MatGetSize(A, NULL, &n));
104 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
105 for (i = rstart; i < rend; ++i) {
106 val = diagfunc(i, n);
107 PetscCall(MatSetValues(A, 1, &i, 1, &i, &val, INSERT_VALUES));
108 }
109 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
110 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
111 PetscFunctionReturn(PETSC_SUCCESS);
112 }
113
main(int argc,char ** argv)114 int main(int argc, char **argv)
115 {
116 PetscInt n = 10000, its, dfid = 1;
117 Vec x, b, u;
118 Mat A;
119 KSP ksp;
120 PC pc, pcnoise;
121 PCNoise_Ctx ctx = {0, NULL};
122 PetscReal eta = 0.1, norm;
123 PetscScalar (*diagfunc)(PetscInt, PetscInt);
124
125 PetscFunctionBeginUser;
126 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
127 /* Process command line options */
128 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
129 PetscCall(PetscOptionsGetReal(NULL, NULL, "-eta", &eta, NULL));
130 PetscCall(PetscOptionsGetInt(NULL, NULL, "-diagfunc", &dfid, NULL));
131 switch (dfid) {
132 case 1:
133 diagfunc = diagFunc1;
134 break;
135 case 2:
136 diagfunc = diagFunc2;
137 break;
138 case 3:
139 diagfunc = diagFunc3;
140 break;
141 default:
142 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unrecognized diagfunc option");
143 }
144
145 /* Create a diagonal matrix with a given distribution of diagonal elements */
146 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
147 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
148 PetscCall(MatSetFromOptions(A));
149 PetscCall(MatSetUp(A));
150 PetscCall(AssembleDiagonalMatrix(A, diagfunc));
151
152 /* Allocate vectors and manufacture an exact solution and rhs */
153 PetscCall(MatCreateVecs(A, &x, NULL));
154 PetscCall(PetscObjectSetName((PetscObject)x, "Computed Solution"));
155 PetscCall(MatCreateVecs(A, &b, NULL));
156 PetscCall(PetscObjectSetName((PetscObject)b, "RHS"));
157 PetscCall(MatCreateVecs(A, &u, NULL));
158 PetscCall(PetscObjectSetName((PetscObject)u, "Reference Solution"));
159 PetscCall(VecSet(u, 1.0));
160 PetscCall(MatMult(A, u, b));
161
162 /* Create a KSP object */
163 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
164 PetscCall(KSPSetOperators(ksp, A, A));
165
166 /* Set up a composite preconditioner */
167 PetscCall(KSPGetPC(ksp, &pc));
168 PetscCall(PCSetType(pc, PCCOMPOSITE)); /* default composite with single Identity PC */
169 PetscCall(PCCompositeSetType(pc, PC_COMPOSITE_ADDITIVE));
170 PetscCall(PCCompositeAddPCType(pc, PCNONE));
171 if (eta > 0) {
172 PetscCall(PCCompositeAddPCType(pc, PCSHELL));
173 PetscCall(PCCompositeGetPC(pc, 1, &pcnoise));
174 ctx.eta = eta;
175 PetscCall(PCShellSetContext(pcnoise, &ctx));
176 PetscCall(PCShellSetApply(pcnoise, PCApply_Noise));
177 PetscCall(PCShellSetSetUp(pcnoise, PCSetup_Noise));
178 PetscCall(PCShellSetDestroy(pcnoise, PCDestroy_Noise));
179 PetscCall(PCShellSetName(pcnoise, "Noise PC"));
180 }
181
182 /* Set KSP from options (this can override the PC just defined) */
183 PetscCall(KSPSetFromOptions(ksp));
184
185 /* Solve */
186 PetscCall(KSPSolve(ksp, b, x));
187
188 /* Compute error */
189 PetscCall(VecAXPY(x, -1.0, u));
190 PetscCall(PetscObjectSetName((PetscObject)x, "Error"));
191 PetscCall(VecNorm(x, NORM_2, &norm));
192 PetscCall(KSPGetIterationNumber(ksp, &its));
193 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
194
195 /* Destroy objects and finalize */
196 PetscCall(KSPDestroy(&ksp));
197 PetscCall(MatDestroy(&A));
198 PetscCall(VecDestroy(&x));
199 PetscCall(VecDestroy(&b));
200 PetscCall(VecDestroy(&u));
201 PetscCall(PetscFinalize());
202 return 0;
203 }
204
205 /*TEST
206
207 build:
208 requires: !complex !single
209
210 test:
211 suffix: 1
212 nsize: 2
213 args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 1 -ksp_type fcg -ksp_fcg_mmax 1 -eta 0.1
214
215 test:
216 suffix: 1_eigs
217 nsize: 2
218 args: -ksp_rtol 1e-6 -diagfunc 1 -ksp_type fcg -ksp_fcg_mmax 1 -eta 0.1 -ksp_monitor_singular_value
219
220 test:
221 suffix: 2
222 nsize: 2
223 args: -ksp_monitor_short -diagfunc 3 -ksp_type fcg -ksp_fcg_mmax 10000 -eta 0.3333
224
225 test:
226 suffix: 3
227 nsize: 3
228 args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 2 -ksp_type fgmres -eta 0.1
229
230 test:
231 suffix: 4
232 nsize: 2
233 args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 1 -ksp_type pipefcg -ksp_pipefcg_mmax 1 -eta 0.1
234
235 test:
236 suffix: 5
237 nsize: 2
238 args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 3 -ksp_type pipefcg -ksp_pipefcg_mmax 10000 -eta 0.1
239
240 test:
241 suffix: 6
242 nsize: 4
243 args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 3 -ksp_type fcg -ksp_fcg_mmax 10000 -eta 0 -pc_type ksp -ksp_ksp_type cg -ksp_pc_type none -ksp_ksp_rtol 1e-1 -ksp_ksp_max_it 5 -ksp_ksp_converged_reason
244
245 test:
246 suffix: 7
247 nsize: 4
248 args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 3 -ksp_type pipefcg -ksp_pipefcg_mmax 10000 -eta 0 -pc_type ksp -ksp_ksp_type cg -ksp_pc_type none -ksp_ksp_rtol 1e-1 -ksp_ksp_max_it 5 -ksp_ksp_converged_reason
249
250 test:
251 suffix: 8
252 nsize: 2
253 args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 1 -ksp_type pipefgmres -pc_type ksp -ksp_ksp_type cg -ksp_pc_type none -ksp_ksp_rtol 1e-2 -ksp_ksp_converged_reason
254
255 test:
256 suffix: 9
257 nsize: 2
258 args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 1 -ksp_type pipefgmres -pc_type ksp -ksp_ksp_type cg -ksp_pc_type none -ksp_ksp_rtol 1e-2 -ksp_ksp_converged_reason
259
260 TEST*/
261