Actual source code: pipecg.c

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

  3: /*
  4:      KSPSetUp_PIPECG - Sets up the workspace needed by the PIPECG 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_PIPECG(KSP ksp)
 10: {
 11:   PetscFunctionBegin;
 12:   /* get work vectors needed by PIPECG */
 13:   PetscCall(KSPSetWorkVecs(ksp, 9));
 14:   PetscFunctionReturn(PETSC_SUCCESS);
 15: }

 17: /*
 18:  KSPSolve_PIPECG - This routine actually applies the pipelined conjugate gradient method
 19: */
 20: static PetscErrorCode KSPSolve_PIPECG(KSP ksp)
 21: {
 22:   PetscInt    i;
 23:   PetscScalar alpha = 0.0, beta = 0.0, gamma = 0.0, gammaold = 0.0, delta = 0.0;
 24:   PetscReal   dp = 0.0;
 25:   Vec         X, B, Z, P, W, Q, U, M, N, R, S;
 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:   R = 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];
 42:   M = ksp->work[7];
 43:   S = ksp->work[8];

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

 47:   ksp->its = 0;
 48:   if (!ksp->guess_zero) {
 49:     PetscCall(KSP_MatMult(ksp, Amat, X, R)); /*     r <- b - Ax     */
 50:     PetscCall(VecAYPX(R, -1.0, B));
 51:   } else {
 52:     PetscCall(VecCopy(B, R)); /*     r <- b (x is 0) */
 53:   }

 55:   PetscCall(KSP_PCApply(ksp, R, U)); /*     u <- Br   */

 57:   switch (ksp->normtype) {
 58:   case KSP_NORM_PRECONDITIONED:
 59:     PetscCall(VecNormBegin(U, NORM_2, &dp)); /*     dp <- u'*u = e'*A'*B'*B*A'*e'     */
 60:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)U)));
 61:     PetscCall(KSP_MatMult(ksp, Amat, U, W)); /*     w <- Au   */
 62:     PetscCall(VecNormEnd(U, NORM_2, &dp));
 63:     break;
 64:   case KSP_NORM_UNPRECONDITIONED:
 65:     PetscCall(VecNormBegin(R, NORM_2, &dp)); /*     dp <- r'*r = e'*A'*A*e            */
 66:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R)));
 67:     PetscCall(KSP_MatMult(ksp, Amat, U, W)); /*     w <- Au   */
 68:     PetscCall(VecNormEnd(R, NORM_2, &dp));
 69:     break;
 70:   case KSP_NORM_NATURAL:
 71:     PetscCall(VecDotBegin(R, U, &gamma)); /*     gamma <- u'*r       */
 72:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R)));
 73:     PetscCall(KSP_MatMult(ksp, Amat, U, W)); /*     w <- Au   */
 74:     PetscCall(VecDotEnd(R, U, &gamma));
 75:     KSPCheckDot(ksp, gamma);
 76:     dp = PetscSqrtReal(PetscAbsScalar(gamma)); /*     dp <- r'*u = r'*B*r = e'*A'*B*A*e */
 77:     break;
 78:   case KSP_NORM_NONE:
 79:     PetscCall(KSP_MatMult(ksp, Amat, U, W));
 80:     dp = 0.0;
 81:     break;
 82:   default:
 83:     SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
 84:   }
 85:   PetscCall(KSPLogResidualHistory(ksp, dp));
 86:   PetscCall(KSPMonitor(ksp, 0, dp));
 87:   ksp->rnorm = dp;
 88:   PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
 89:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

 91:   i = 0;
 92:   do {
 93:     if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 94:       PetscCall(VecNormBegin(R, NORM_2, &dp));
 95:     } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
 96:       PetscCall(VecNormBegin(U, NORM_2, &dp));
 97:     }
 98:     if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) PetscCall(VecDotBegin(R, U, &gamma));
 99:     PetscCall(VecDotBegin(W, U, &delta));
100:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R)));

102:     PetscCall(KSP_PCApply(ksp, W, M));       /*   m <- Bw       */
103:     PetscCall(KSP_MatMult(ksp, Amat, M, N)); /*   n <- Am       */

105:     if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
106:       PetscCall(VecNormEnd(R, NORM_2, &dp));
107:     } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
108:       PetscCall(VecNormEnd(U, NORM_2, &dp));
109:     }
110:     if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) PetscCall(VecDotEnd(R, U, &gamma));
111:     PetscCall(VecDotEnd(W, U, &delta));

113:     if (i > 0) {
114:       if (ksp->normtype == KSP_NORM_NATURAL) dp = PetscSqrtReal(PetscAbsScalar(gamma));
115:       else if (ksp->normtype == KSP_NORM_NONE) dp = 0.0;

117:       ksp->rnorm = dp;
118:       PetscCall(KSPLogResidualHistory(ksp, dp));
119:       PetscCall(KSPMonitor(ksp, i, dp));
120:       PetscCall((*ksp->converged)(ksp, i, dp, &ksp->reason, ksp->cnvP));
121:       if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
122:     }

124:     if (i == 0) {
125:       alpha = gamma / delta;
126:       PetscCall(VecCopy(N, Z)); /*     z <- n          */
127:       PetscCall(VecCopy(M, Q)); /*     q <- m          */
128:       PetscCall(VecCopy(U, P)); /*     p <- u          */
129:       PetscCall(VecCopy(W, S)); /*     s <- w          */
130:     } else {
131:       beta  = gamma / gammaold;
132:       alpha = gamma / (delta - beta / alpha * gamma);
133:       PetscCall(VecAYPX(Z, beta, N)); /*     z <- n + beta * z   */
134:       PetscCall(VecAYPX(Q, beta, M)); /*     q <- m + beta * q   */
135:       PetscCall(VecAYPX(P, beta, U)); /*     p <- u + beta * p   */
136:       PetscCall(VecAYPX(S, beta, W)); /*     s <- w + beta * s   */
137:     }
138:     PetscCall(VecAXPY(X, alpha, P));  /*     x <- x + alpha * p   */
139:     PetscCall(VecAXPY(U, -alpha, Q)); /*     u <- u - alpha * q   */
140:     PetscCall(VecAXPY(W, -alpha, Z)); /*     w <- w - alpha * z   */
141:     PetscCall(VecAXPY(R, -alpha, S)); /*     r <- r - alpha * s   */
142:     gammaold = gamma;
143:     i++;
144:     ksp->its = i;

146:     /* if (i%50 == 0) { */
147:     /*   PetscCall(KSP_MatMult(ksp,Amat,X,R));            /\*     w <- b - Ax     *\/ */
148:     /*   PetscCall(VecAYPX(R,-1.0,B)); */
149:     /*   PetscCall(KSP_PCApply(ksp,R,U)); */
150:     /*   PetscCall(KSP_MatMult(ksp,Amat,U,W)); */
151:     /* } */

153:   } while (i <= ksp->max_it);
154:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }

158: PETSC_INTERN PetscErrorCode KSPBuildResidual_CG(KSP, Vec, Vec, Vec *);

160: /*MC
161:    KSPPIPECG - Pipelined conjugate gradient method {cite}`ghyselsvanroose2014`. [](sec_pipelineksp)

163:    Level: intermediate

165:    Notes:
166:    This method has only a single non-blocking reduction per iteration, compared to 2 blocking for standard `KSPPCG`.  The
167:    non-blocking reduction is overlapped by the matrix-vector product and preconditioner application.

169:    See also `KSPPIPECR`, where the reduction is only overlapped with the matrix-vector product and `KSPGROPPCG`

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

174:    Contributed by:
175:    Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders

177: .seealso: [](ch_ksp), [](doc_faq_pipelined), [](sec_pipelineksp), `KSPCreate()`, `KSPSetType()`, `KSPPIPECG2`, `KSPPIPECR`, `KSPGROPPCG`, `KSPPGMRES`, `KSPCG`, `KSPCGUseSingleReduction()`
178: M*/
179: PETSC_EXTERN PetscErrorCode KSPCreate_PIPECG(KSP ksp)
180: {
181:   PetscFunctionBegin;
182:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
183:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
184:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
185:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

187:   ksp->ops->setup          = KSPSetUp_PIPECG;
188:   ksp->ops->solve          = KSPSolve_PIPECG;
189:   ksp->ops->destroy        = KSPDestroyDefault;
190:   ksp->ops->view           = NULL;
191:   ksp->ops->setfromoptions = NULL;
192:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
193:   ksp->ops->buildresidual  = KSPBuildResidual_CG;
194:   PetscFunctionReturn(PETSC_SUCCESS);
195: }