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