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;
43: 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: VecSetRandom(xout, ctx->random);
48: VecNorm(xin, NORM_2, &nrmin);
49: VecNorm(xout, NORM_2, &nrmnoise);
50: VecScale(xout, ctx->eta * (nrmin / nrmnoise));
51: return 0;
52: }
54: PetscErrorCode PCSetup_Noise(PC pc)
55: {
56: PCNoise_Ctx *ctx;
59: PCShellGetContext(pc, &ctx);
60: PetscRandomCreate(PETSC_COMM_WORLD, &ctx->random);
61: PetscRandomSetInterval(ctx->random, -1.0, 1.0);
62: return 0;
63: }
65: PetscErrorCode PCDestroy_Noise(PC pc)
66: {
67: PCNoise_Ctx *ctx;
70: PCShellGetContext(pc, &ctx);
71: PetscRandomDestroy(&ctx->random);
72: return 0;
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;
103: MatGetSize(A, NULL, &n);
104: MatGetOwnershipRange(A, &rstart, &rend);
105: for (i = rstart; i < rend; ++i) {
106: val = diagfunc(i, n);
107: MatSetValues(A, 1, &i, 1, &i, &val, INSERT_VALUES);
108: }
109: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
110: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
111: return 0;
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);
126: PetscInitialize(&argc, &argv, (char *)0, help);
127: /* Process command line options */
128: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
129: PetscOptionsGetReal(NULL, NULL, "-eta", &eta, NULL);
130: 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: MatCreate(PETSC_COMM_WORLD, &A);
147: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
148: MatSetFromOptions(A);
149: MatSetUp(A);
150: AssembleDiagonalMatrix(A, diagfunc);
152: /* Allocate vectors and manufacture an exact solution and rhs */
153: MatCreateVecs(A, &x, NULL);
154: PetscObjectSetName((PetscObject)x, "Computed Solution");
155: MatCreateVecs(A, &b, NULL);
156: PetscObjectSetName((PetscObject)b, "RHS");
157: MatCreateVecs(A, &u, NULL);
158: PetscObjectSetName((PetscObject)u, "Reference Solution");
159: VecSet(u, 1.0);
160: MatMult(A, u, b);
162: /* Create a KSP object */
163: KSPCreate(PETSC_COMM_WORLD, &ksp);
164: KSPSetOperators(ksp, A, A);
166: /* Set up a composite preconditioner */
167: KSPGetPC(ksp, &pc);
168: PCSetType(pc, PCCOMPOSITE); /* default composite with single Identity PC */
169: PCCompositeSetType(pc, PC_COMPOSITE_ADDITIVE);
170: PCCompositeAddPCType(pc, PCNONE);
171: if (eta > 0) {
172: PCCompositeAddPCType(pc, PCSHELL);
173: PCCompositeGetPC(pc, 1, &pcnoise);
174: ctx.eta = eta;
175: PCShellSetContext(pcnoise, &ctx);
176: PCShellSetApply(pcnoise, PCApply_Noise);
177: PCShellSetSetUp(pcnoise, PCSetup_Noise);
178: PCShellSetDestroy(pcnoise, PCDestroy_Noise);
179: PCShellSetName(pcnoise, "Noise PC");
180: }
182: /* Set KSP from options (this can override the PC just defined) */
183: KSPSetFromOptions(ksp);
185: /* Solve */
186: KSPSolve(ksp, b, x);
188: /* Compute error */
189: VecAXPY(x, -1.0, u);
190: PetscObjectSetName((PetscObject)x, "Error");
191: VecNorm(x, NORM_2, &norm);
192: KSPGetIterationNumber(ksp, &its);
193: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
195: /* Destroy objects and finalize */
196: KSPDestroy(&ksp);
197: MatDestroy(&A);
198: VecDestroy(&x);
199: VecDestroy(&b);
200: VecDestroy(&u);
201: PetscFinalize();
202: return 0;
203: }
205: /*TEST
207: build:
208: requires: !complex !single
210: test:
211: nsize: 2
212: args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 1 -ksp_type fcg -ksp_fcg_mmax 1 -eta 0.1
214: test:
215: suffix: 2
216: nsize: 2
217: args: -ksp_monitor_short -diagfunc 3 -ksp_type fcg -ksp_fcg_mmax 10000 -eta 0.3333
219: test:
220: suffix: 3
221: nsize: 3
222: args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 2 -ksp_type fgmres -eta 0.1
224: test:
225: suffix: 4
226: nsize: 2
227: args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 1 -ksp_type pipefcg -ksp_pipefcg_mmax 1 -eta 0.1
229: test:
230: suffix: 5
231: nsize: 2
232: args: -ksp_monitor_short -ksp_rtol 1e-6 -diagfunc 3 -ksp_type pipefcg -ksp_pipefcg_mmax 10000 -eta 0.1
234: test:
235: suffix: 6
236: nsize: 4
237: 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
239: test:
240: suffix: 7
241: nsize: 4
242: 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
244: test:
245: suffix: 8
246: nsize: 2
247: 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
249: test:
250: suffix: 9
251: nsize: 2
252: 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: TEST*/