Actual source code: groppcg.c

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

  3: /*
  4:  KSPSetUp_GROPPCG - Sets up the workspace needed by the GROPPCG 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_GROPPCG(KSP ksp)
 10: {
 11:   PetscFunctionBegin;
 12:   PetscCall(KSPSetWorkVecs(ksp, 6));
 13:   PetscFunctionReturn(PETSC_SUCCESS);
 14: }

 16: /*
 17:  KSPSolve_GROPPCG

 19:  Input Parameter:
 20:  .     ksp - the Krylov space object that was set to use conjugate gradient, by, for
 21:              example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG);
 22: */
 23: static PetscErrorCode KSPSolve_GROPPCG(KSP ksp)
 24: {
 25:   PetscInt    i;
 26:   PetscScalar alpha, beta = 0.0, gamma, gammaNew, t;
 27:   PetscReal   dp = 0.0;
 28:   Vec         x, b, r, p, s, S, z, Z;
 29:   Mat         Amat, Pmat;
 30:   PetscBool   diagonalscale;

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

 36:   x = ksp->vec_sol;
 37:   b = ksp->vec_rhs;
 38:   r = ksp->work[0];
 39:   p = ksp->work[1];
 40:   s = ksp->work[2];
 41:   S = ksp->work[3];
 42:   z = ksp->work[4];
 43:   Z = ksp->work[5];

 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, z));    /*     z <- Br   */
 56:   PetscCall(VecCopy(z, p));             /*     p <- z    */
 57:   PetscCall(VecDotBegin(r, z, &gamma)); /*     gamma <- z'*r       */
 58:   PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r)));
 59:   PetscCall(KSP_MatMult(ksp, Amat, p, s)); /*     s <- Ap   */
 60:   PetscCall(VecDotEnd(r, z, &gamma));      /*     gamma <- z'*r       */

 62:   switch (ksp->normtype) {
 63:   case KSP_NORM_PRECONDITIONED:
 64:     /* This could be merged with the computation of gamma above */
 65:     PetscCall(VecNorm(z, NORM_2, &dp)); /*     dp <- z'*z = e'*A'*B'*B*A'*e'     */
 66:     break;
 67:   case KSP_NORM_UNPRECONDITIONED:
 68:     /* This could be merged with the computation of gamma above */
 69:     PetscCall(VecNorm(r, NORM_2, &dp)); /*     dp <- r'*r = e'*A'*A*e            */
 70:     break;
 71:   case KSP_NORM_NATURAL:
 72:     KSPCheckDot(ksp, gamma);
 73:     dp = PetscSqrtReal(PetscAbsScalar(gamma)); /*     dp <- r'*z = r'*B*r = e'*A'*B*A*e */
 74:     break;
 75:   case KSP_NORM_NONE:
 76:     dp = 0.0;
 77:     break;
 78:   default:
 79:     SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
 80:   }
 81:   PetscCall(KSPLogResidualHistory(ksp, dp));
 82:   PetscCall(KSPMonitor(ksp, 0, dp));
 83:   ksp->rnorm = dp;
 84:   PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
 85:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

 87:   i = 0;
 88:   do {
 89:     ksp->its = i + 1;
 90:     i++;

 92:     PetscCall(VecDotBegin(p, s, &t));
 93:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)p)));

 95:     PetscCall(KSP_PCApply(ksp, s, S)); /*   S <- Bs       */

 97:     PetscCall(VecDotEnd(p, s, &t));

 99:     alpha = gamma / t;
100:     PetscCall(VecAXPY(x, alpha, p));  /*     x <- x + alpha * p   */
101:     PetscCall(VecAXPY(r, -alpha, s)); /*     r <- r - alpha * s   */
102:     PetscCall(VecAXPY(z, -alpha, S)); /*     z <- z - alpha * S   */

104:     if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
105:       PetscCall(VecNormBegin(r, NORM_2, &dp));
106:     } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
107:       PetscCall(VecNormBegin(z, NORM_2, &dp));
108:     }
109:     PetscCall(VecDotBegin(r, z, &gammaNew));
110:     PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r)));

112:     PetscCall(KSP_MatMult(ksp, Amat, z, Z)); /*   Z <- Az       */

114:     if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
115:       PetscCall(VecNormEnd(r, NORM_2, &dp));
116:     } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
117:       PetscCall(VecNormEnd(z, NORM_2, &dp));
118:     }
119:     PetscCall(VecDotEnd(r, z, &gammaNew));

121:     if (ksp->normtype == KSP_NORM_NATURAL) {
122:       KSPCheckDot(ksp, gammaNew);
123:       dp = PetscSqrtReal(PetscAbsScalar(gammaNew)); /*     dp <- r'*z = r'*B*r = e'*A'*B*A*e */
124:     } else if (ksp->normtype == KSP_NORM_NONE) {
125:       dp = 0.0;
126:     }
127:     ksp->rnorm = dp;
128:     PetscCall(KSPLogResidualHistory(ksp, dp));
129:     PetscCall(KSPMonitor(ksp, i, dp));
130:     PetscCall((*ksp->converged)(ksp, i, dp, &ksp->reason, ksp->cnvP));
131:     if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

133:     beta  = gammaNew / gamma;
134:     gamma = gammaNew;
135:     PetscCall(VecAYPX(p, beta, z)); /*     p <- z + beta * p   */
136:     PetscCall(VecAYPX(s, beta, Z)); /*     s <- Z + beta * s   */

138:   } while (i < ksp->max_it);

140:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
141:   PetscFunctionReturn(PETSC_SUCCESS);
142: }

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

146: /*MC
147:    KSPGROPPCG - A pipelined conjugate gradient method developed by Bill Gropp {cite}`eller2016scalable`. [](sec_pipelineksp)

149:    Level: intermediate

151:    Notes:
152:    This method has two reductions, one of which is overlapped with the matrix-vector product and one of which is
153:    overlapped with the preconditioner.

155:    See also `KSPPIPECG`, which has only a single reduction that overlaps both the matrix-vector product and the preconditioner.

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

160:    Contributed by:
161:    Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders

163: .seealso: [](ch_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPPIPECG2()`, `KSPSetType()`, `KSPPIPECG`, `KSPPIPECR`, `KSPPGMRES`, `KSPCG`, `KSPCGUseSingleReduction()`
164: M*/

166: PETSC_EXTERN PetscErrorCode KSPCreate_GROPPCG(KSP ksp)
167: {
168:   PetscFunctionBegin;
169:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
170:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
171:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
172:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

174:   ksp->ops->setup          = KSPSetUp_GROPPCG;
175:   ksp->ops->solve          = KSPSolve_GROPPCG;
176:   ksp->ops->destroy        = KSPDestroyDefault;
177:   ksp->ops->view           = NULL;
178:   ksp->ops->setfromoptions = NULL;
179:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
180:   ksp->ops->buildresidual  = KSPBuildResidual_CG;
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }