Actual source code: pipecr.c

      1: #include <petsc/private/kspimpl.h>

      3: /*
      4:      KSPSetUp_PIPECR - Sets up the workspace needed by the PIPECR method.

      6:       This is called once, usually automatically by KSPSolve() or KSPSetUp()
      7:      but can be called directly by KSPSetUp()
      8: */
      9: static PetscErrorCode KSPSetUp_PIPECR(KSP ksp)
     10: {
     11:   PetscFunctionBegin;
     12:   /* get work vectors needed by PIPECR */
     13:   PetscCall(KSPSetWorkVecs(ksp, 7));
     14:   PetscFunctionReturn(PETSC_SUCCESS);
     15: }

     17: /*
     18:  KSPSolve_PIPECR - This routine actually applies the pipelined conjugate residual method
     19: */
     20: static PetscErrorCode KSPSolve_PIPECR(KSP ksp)
     21: {
     22:   PetscInt    i;
     23:   PetscScalar alpha = 0.0, beta = 0.0, gamma, gammaold = 0.0, delta;
     24:   PetscReal   dp = 0.0;
     25:   Vec         X, B, Z, P, W, Q, U, M, N;
     26:   Mat         Amat, Pmat;
     27:   PetscBool   diagonalscale;

     29:   PetscFunctionBegin;
     30:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
     31:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

     33:   X = ksp->vec_sol;
     34:   B = ksp->vec_rhs;
     35:   M = ksp->work[0];
     36:   Z = ksp->work[1];
     37:   P = ksp->work[2];
     38:   N = ksp->work[3];
     39:   W = ksp->work[4];
     40:   Q = ksp->work[5];
     41:   U = ksp->work[6];

     43:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));

     45:   ksp->its = 0;
     46:   /* we don't have an R vector, so put the (unpreconditioned) residual in w for now */
     47:   if (!ksp->guess_zero) {
     48:     PetscCall(KSP_MatMult(ksp, Amat, X, W)); /*     w <- b - Ax     */
     49:     PetscCall(VecAYPX(W, -1.0, B));
     50:   } else {
     51:     PetscCall(VecCopy(B, W)); /*     w <- b (x is 0) */
     52:   }
     53:   PetscCall(KSP_PCApply(ksp, W, U)); /*     u <- Bw   */

     55:   switch (ksp->normtype) {
     56:   case KSP_NORM_PRECONDITIONED:
     57:     PetscCall(VecNormBegin(U, NORM_2, &dp)); /*     dp <- u'*u = e'*A'*B'*B*A'*e'     */
     58:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)U)));
     59:     PetscCall(KSP_MatMult(ksp, Amat, U, W)); /*     w <- Au   */
     60:     PetscCall(VecNormEnd(U, NORM_2, &dp));
     61:     break;
     62:   case KSP_NORM_NONE:
     63:     PetscCall(KSP_MatMult(ksp, Amat, U, W));
     64:     dp = 0.0;
     65:     break;
     66:   default:
     67:     SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
     68:   }
     69:   PetscCall(KSPLogResidualHistory(ksp, dp));
     70:   PetscCall(KSPMonitor(ksp, 0, dp));
     71:   ksp->rnorm = dp;
     72:   PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
     73:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

     75:   i = 0;
     76:   do {
     77:     PetscCall(KSP_PCApply(ksp, W, M)); /*   m <- Bw       */

     79:     if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) PetscCall(VecNormBegin(U, NORM_2, &dp));
     80:     PetscCall(VecDotBegin(W, U, &gamma));
     81:     PetscCall(VecDotBegin(M, W, &delta));
     82:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)U)));

     84:     PetscCall(KSP_MatMult(ksp, Amat, M, N)); /*   n <- Am       */

     86:     if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) PetscCall(VecNormEnd(U, NORM_2, &dp));
     87:     PetscCall(VecDotEnd(W, U, &gamma));
     88:     PetscCall(VecDotEnd(M, W, &delta));

     90:     if (i > 0) {
     91:       if (ksp->normtype == KSP_NORM_NONE) dp = 0.0;
     92:       ksp->rnorm = dp;
     93:       PetscCall(KSPLogResidualHistory(ksp, dp));
     94:       PetscCall(KSPMonitor(ksp, i, dp));
     95:       PetscCall((*ksp->converged)(ksp, i, dp, &ksp->reason, ksp->cnvP));
     96:       if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
     97:     }

     99:     if (i == 0) {
    100:       alpha = gamma / delta;
    101:       PetscCall(VecCopy(N, Z)); /*     z <- n          */
    102:       PetscCall(VecCopy(M, Q)); /*     q <- m          */
    103:       PetscCall(VecCopy(U, P)); /*     p <- u          */
    104:     } else {
    105:       beta  = gamma / gammaold;
    106:       alpha = gamma / (delta - beta / alpha * gamma);
    107:       PetscCall(VecAYPX(Z, beta, N)); /*     z <- n + beta * z   */
    108:       PetscCall(VecAYPX(Q, beta, M)); /*     q <- m + beta * q   */
    109:       PetscCall(VecAYPX(P, beta, U)); /*     p <- u + beta * p   */
    110:     }
    111:     PetscCall(VecAXPY(X, alpha, P));  /*     x <- x + alpha * p   */
    112:     PetscCall(VecAXPY(U, -alpha, Q)); /*     u <- u - alpha * q   */
    113:     PetscCall(VecAXPY(W, -alpha, Z)); /*     w <- w - alpha * z   */
    114:     gammaold = gamma;
    115:     i++;
    116:     ksp->its = i;

    118:     /* if (i%50 == 0) { */
    119:     /*   PetscCall(KSP_MatMult(ksp,Amat,X,W));            /\*     w <- b - Ax     *\/ */
    120:     /*   PetscCall(VecAYPX(W,-1.0,B)); */
    121:     /*   PetscCall(KSP_PCApply(ksp,W,U)); */
    122:     /*   PetscCall(KSP_MatMult(ksp,Amat,U,W)); */
    123:     /* } */

    125:   } while (i <= ksp->max_it);
    126:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    127:   PetscFunctionReturn(PETSC_SUCCESS);
    128: }

    130: /*MC
    131:    KSPPIPECR - Pipelined conjugate residual method {cite}`ghyselsvanroose2014`. [](sec_pipelineksp)

    133:    Level: intermediate

    135:    Notes:
    136:    This method has only a single non-blocking reduction per iteration, compared to 2 for standard `KSPCR`.  The
    137:    non-blocking reduction is overlapped by the matrix-vector product, but not the preconditioner application.

    139:    See also `KSPPIPECG`, where the reduction is overlapped with the matrix-vector product.

    141:    MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
    142:    See [](doc_faq_pipelined)

    144:    Contributed by:
    145:    Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders

    147: .seealso: [](ch_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPSetType()`, `KSPPIPECG`, `KSPGROPPCG`, `KSPPGMRES`, `KSPCG`, `KSPCGUseSingleReduction()`
    148: M*/

    150: PETSC_EXTERN PetscErrorCode KSPCreate_PIPECR(KSP ksp)
    151: {
    152:   PetscFunctionBegin;
    153:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
    154:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

    156:   ksp->ops->setup          = KSPSetUp_PIPECR;
    157:   ksp->ops->solve          = KSPSolve_PIPECR;
    158:   ksp->ops->destroy        = KSPDestroyDefault;
    159:   ksp->ops->view           = NULL;
    160:   ksp->ops->setfromoptions = NULL;
    161:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    162:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    163:   PetscFunctionReturn(PETSC_SUCCESS);
    164: }