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*/