Actual source code: ex60.c
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";
8: /*
9: Solve (in parallel) a diagonal linear system.
11: Using PCShell, we define a preconditioner which simply adds noise to the residual.
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.
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.
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)
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
26: Contributed by Patrick Sanan
27: */
29: #include <petscksp.h>
31: /* Context to use with our noise PC */
32: typedef struct {
33: PetscReal eta;
34: PetscRandom random;
35: } PCNoise_Ctx;
37: PetscErrorCode PCApply_Noise(PC pc, Vec xin, Vec xout)
38: {
39: PCNoise_Ctx *ctx;
40: PetscReal nrmin, nrmnoise;
42: PetscFunctionBeginUser;
43: PetscCall(PCShellGetContext(pc, &ctx));
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: }
54: PetscErrorCode PCSetup_Noise(PC pc)
55: {
56: PCNoise_Ctx *ctx;
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: }
65: PetscErrorCode PCDestroy_Noise(PC pc)
66: {
67: PCNoise_Ctx *ctx;
69: PetscFunctionBeginUser;
70: PetscCall(PCShellGetContext(pc, &ctx));
71: PetscCall(PetscRandomDestroy(&ctx->random));
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
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: }
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: }
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: }
97: static PetscErrorCode AssembleDiagonalMatrix(Mat A, PetscScalar (*diagfunc)(PetscInt, PetscInt))
98: {
99: PetscInt i, rstart, rend, n;
100: PetscScalar val;
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: }
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);
125: PetscFunctionBeginUser;
126: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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: }
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));
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));
162: /* Create a KSP object */
163: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
164: PetscCall(KSPSetOperators(ksp, A, A));
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: }
182: /* Set KSP from options (this can override the PC just defined) */
183: PetscCall(KSPSetFromOptions(ksp));
185: /* Solve */
186: PetscCall(KSPSolve(ksp, b, x));
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));
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: }
205: /*TEST
207: build:
208: requires: !complex !single
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
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
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
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
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
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
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
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
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
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
260: TEST*/