Actual source code: pipecg.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

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

  4: /*
  5:      KSPSetUp_PIPECG - Sets up the workspace needed by the PIPECG method.

  7:       This is called once, usually automatically by KSPSolve() or KSPSetUp()
  8:      but can be called directly by KSPSetUp()
  9: */
 10: static PetscErrorCode KSPSetUp_PIPECG(KSP ksp)
 11: {

 15:   /* get work vectors needed by PIPECG */
 16:   KSPSetWorkVecs(ksp,9);
 17:   return(0);
 18: }

 20: /*
 21:  KSPSolve_PIPECG - This routine actually applies the pipelined conjugate gradient method
 22: */
 23: static PetscErrorCode  KSPSolve_PIPECG(KSP ksp)
 24: {
 26:   PetscInt       i;
 27:   PetscScalar    alpha = 0.0,beta = 0.0,gamma = 0.0,gammaold = 0.0,delta = 0.0;
 28:   PetscReal      dp    = 0.0;
 29:   Vec            X,B,Z,P,W,Q,U,M,N,R,S;
 30:   Mat            Amat,Pmat;
 31:   PetscBool      diagonalscale;

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

 37:   X = ksp->vec_sol;
 38:   B = ksp->vec_rhs;
 39:   R = ksp->work[0];
 40:   Z = ksp->work[1];
 41:   P = ksp->work[2];
 42:   N = ksp->work[3];
 43:   W = ksp->work[4];
 44:   Q = ksp->work[5];
 45:   U = ksp->work[6];
 46:   M = ksp->work[7];
 47:   S = ksp->work[8];

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

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

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

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

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

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

110:     if (i > 0 && ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
111:       VecNormEnd(R,NORM_2,&dp);
112:     } else if (i > 0 && ksp->normtype == KSP_NORM_PRECONDITIONED) {
113:       VecNormEnd(U,NORM_2,&dp);
114:     }
115:     if (!(i == 0 && ksp->normtype == KSP_NORM_NATURAL)) {
116:       VecDotEnd(R,U,&gamma);
117:     }
118:     VecDotEnd(W,U,&delta);

120:     if (i > 0) {
121:       if (ksp->normtype == KSP_NORM_NATURAL) dp = PetscSqrtReal(PetscAbsScalar(gamma));
122:       else if (ksp->normtype == KSP_NORM_NONE) dp = 0.0;

124:       ksp->rnorm = dp;
125:       KSPLogResidualHistory(ksp,dp);
126:       KSPMonitor(ksp,i,dp);
127:       (*ksp->converged)(ksp,i,dp,&ksp->reason,ksp->cnvP);
128:       if (ksp->reason) return(0);
129:     }

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

153:     /* if (i%50 == 0) { */
154:     /*   KSP_MatMult(ksp,Amat,X,R);            /\*     w <- b - Ax     *\/ */
155:     /*   VecAYPX(R,-1.0,B); */
156:     /*   KSP_PCApply(ksp,R,U); */
157:     /*   KSP_MatMult(ksp,Amat,U,W); */
158:     /* } */

160:   } while (i<=ksp->max_it);
161:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
162:   return(0);
163: }

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

167: /*MC
168:    KSPPIPECG - Pipelined conjugate gradient method.

170:    This method has only a single non-blocking reduction per iteration, compared to 2 blocking for standard CG.  The
171:    non-blocking reduction is overlapped by the matrix-vector product and preconditioner application.

173:    See also KSPPIPECR, where the reduction is only overlapped with the matrix-vector product.

175:    Level: intermediate

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

181:    Contributed by:
182:    Pieter Ghysels, Universiteit Antwerpen, Intel Exascience lab Flanders

184:    Reference:
185:    P. Ghysels and W. Vanroose, "Hiding global synchronization latency in the preconditioned Conjugate Gradient algorithm",
186:    Submitted to Parallel Computing, 2012.

188: .seealso: KSPCreate(), KSPSetType(), KSPPIPECR, KSPGROPPCG, KSPPGMRES, KSPCG, KSPCGUseSingleReduction()
189: M*/
190: PETSC_EXTERN PetscErrorCode KSPCreate_PIPECG(KSP ksp)
191: {

195:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
196:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
197:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
198:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

200:   ksp->ops->setup          = KSPSetUp_PIPECG;
201:   ksp->ops->solve          = KSPSolve_PIPECG;
202:   ksp->ops->destroy        = KSPDestroyDefault;
203:   ksp->ops->view           = NULL;
204:   ksp->ops->setfromoptions = NULL;
205:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
206:   ksp->ops->buildresidual  = KSPBuildResidual_CG;
207:   return(0);
208: }