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:   KSPSetWorkVecs(ksp, 6);
 12:   return 0;
 13: }

 15: /*
 16:  KSPSolve_GROPPCG

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

 31:   PCGetDiagonalScale(ksp->pc, &diagonalscale);

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

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

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

 53:   KSP_PCApply(ksp, r, z);    /*     z <- Br   */
 54:   VecCopy(z, p);             /*     p <- z    */
 55:   VecDotBegin(r, z, &gamma); /*     gamma <- z'*r       */
 56:   PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));
 57:   KSP_MatMult(ksp, Amat, p, s); /*     s <- Ap   */
 58:   VecDotEnd(r, z, &gamma);      /*     gamma <- z'*r       */

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

 85:   i = 0;
 86:   do {
 87:     ksp->its = i + 1;
 88:     i++;

 90:     VecDotBegin(p, s, &t);
 91:     PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)p));

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

 95:     VecDotEnd(p, s, &t);

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

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

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

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

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

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

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

138:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
139:   return 0;
140: }

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

144: /*MC
145:    KSPGROPPCG - A pipelined conjugate gradient method developed by Bill Gropp. [](sec_pipelineksp)

147:    Level: intermediate

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

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

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

158:    Contributed by:
159:    Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders

161:    Reference:
162:    http://www.cs.uiuc.edu/~wgropp/bib/talks/tdata/2012/icerm.pdf

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

167: PETSC_EXTERN PetscErrorCode KSPCreate_GROPPCG(KSP ksp)
168: {
169:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);
170:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
171:   KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
172:   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:   return 0;
182: }