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, 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:   }

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