Actual source code: groppcg.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  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: */
 11: PetscErrorCode KSPSetUp_GROPPCG(KSP ksp)
 12: {

 16:   KSPSetWorkVecs(ksp,6);
 17:   return(0);
 18: }

 20: /*
 21:  KSPSolve_GROPPCG

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

 40:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
 41:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

 43:   x = ksp->vec_sol;
 44:   b = ksp->vec_rhs;
 45:   r = ksp->work[0];
 46:   p = ksp->work[1];
 47:   s = ksp->work[2];
 48:   S = ksp->work[3];
 49:   z = ksp->work[4];
 50:   Z = ksp->work[5];

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

 54:   ksp->its = 0;
 55:   if (!ksp->guess_zero) {
 56:     KSP_MatMult(ksp,Amat,x,r);            /*     r <- b - Ax     */
 57:     VecAYPX(r,-1.0,b);
 58:   } else {
 59:     VecCopy(b,r);                         /*     r <- b (x is 0) */
 60:   }

 62:   KSP_PCApply(ksp,r,z);                   /*     z <- Br   */
 63:   VecCopy(z,p);                           /*     p <- z    */
 64:   VecDotBegin(r,z,&gamma);                  /*     gamma <- z'*r       */
 65:   PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));
 66:   KSP_MatMult(ksp,Amat,p,s);              /*     s <- Ap   */
 67:   VecDotEnd(r,z,&gamma);                  /*     gamma <- z'*r       */

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

 93:   i = 0;
 94:   do {
 95:     ksp->its = i+1;
 96:     i++;

 98:     VecDotBegin(p,s,&t);
 99:     PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)p));

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

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

105:     alpha = gamma / t;
106:     VecAXPY(x, alpha,p);   /*     x <- x + alpha * p   */
107:     VecAXPY(r,-alpha,s);   /*     r <- r - alpha * s   */
108:     VecAXPY(z,-alpha,S);   /*     z <- z - alpha * S   */

110:     if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
111:       VecNormBegin(r,NORM_2,&dp);
112:     } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
113:       VecNormBegin(z,NORM_2,&dp);
114:     }
115:     VecDotBegin(r,z,&gammaNew);
116:     PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)r));

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

120:     if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
121:       VecNormEnd(r,NORM_2,&dp);
122:     } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
123:       VecNormEnd(z,NORM_2,&dp);
124:     }
125:     VecDotEnd(r,z,&gammaNew);

127:     if (ksp->normtype == KSP_NORM_NATURAL) {
128:       KSPCheckDot(ksp,gammaNew);
129:       dp = PetscSqrtReal(PetscAbsScalar(gammaNew));                  /*     dp <- r'*z = r'*B*r = e'*A'*B*A*e */
130:     } else if (ksp->normtype == KSP_NORM_NONE) {
131:       dp = 0.0;
132:     }
133:     ksp->rnorm = dp;
134:     KSPLogResidualHistory(ksp,dp);
135:     KSPMonitor(ksp,i,dp);
136:     (*ksp->converged)(ksp,i,dp,&ksp->reason,ksp->cnvP);
137:     if (ksp->reason) break;

139:     beta  = gammaNew / gamma;
140:     gamma = gammaNew;
141:     VecAYPX(p,beta,z);   /*     p <- z + beta * p   */
142:     VecAYPX(s,beta,Z);   /*     s <- Z + beta * s   */

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

146:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
147:   return(0);
148: }

150: /*MC
151:    KSPGROPPCG - A pipelined conjugate gradient method from Bill Gropp

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

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

158:    Level: intermediate

160:    Notes:
161:    MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
162:    See the FAQ on the PETSc website for details.

164:    Contributed by:
165:    Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders

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

170: .seealso: KSPCreate(), KSPSetType(), KSPPIPECG, KSPPIPECR, KSPPGMRES, KSPCG, KSPCGUseSingleReduction()
171: M*/

175: PETSC_EXTERN PetscErrorCode KSPCreate_GROPPCG(KSP ksp)
176: {

180:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
181:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
182:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

184:   ksp->ops->setup          = KSPSetUp_GROPPCG;
185:   ksp->ops->solve          = KSPSolve_GROPPCG;
186:   ksp->ops->destroy        = KSPDestroyDefault;
187:   ksp->ops->view           = 0;
188:   ksp->ops->setfromoptions = 0;
189:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
190:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
191:   return(0);
192: }