Actual source code: cr.c

petsc-master 2017-03-28
Report Typos and Errors

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

  4: static PetscErrorCode KSPSetUp_CR(KSP ksp)
  5: {

  9:   if (ksp->pc_side == PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"no right preconditioning for KSPCR");
 10:   else if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"no symmetric preconditioning for KSPCR");
 11:   KSPSetWorkVecs(ksp,6);
 12:   return(0);
 13: }

 15: static PetscErrorCode  KSPSolve_CR(KSP ksp)
 16: {
 18:   PetscInt       i = 0;
 19:   PetscReal      dp;
 20:   PetscScalar    ai, bi;
 21:   PetscScalar    apq,btop, bbot;
 22:   Vec            X,B,R,RT,P,AP,ART,Q;
 23:   Mat            Amat, Pmat;

 26:   X   = ksp->vec_sol;
 27:   B   = ksp->vec_rhs;
 28:   R   = ksp->work[0];
 29:   RT  = ksp->work[1];
 30:   P   = ksp->work[2];
 31:   AP  = ksp->work[3];
 32:   ART = ksp->work[4];
 33:   Q   = ksp->work[5];

 35:   /* R is the true residual norm, RT is the preconditioned residual norm */
 36:   PCGetOperators(ksp->pc,&Amat,&Pmat);
 37:   if (!ksp->guess_zero) {
 38:     KSP_MatMult(ksp,Amat,X,R);     /*   R <- A*X           */
 39:     VecAYPX(R,-1.0,B);            /*   R <- B-R == B-A*X  */
 40:   } else {
 41:     VecCopy(B,R);                  /*   R <- B (X is 0)    */
 42:   }
 43:   KSP_PCApply(ksp,R,P);     /*   P   <- B*R         */
 44:   KSP_MatMult(ksp,Amat,P,AP);      /*   AP  <- A*P         */
 45:   VecCopy(P,RT);                   /*   RT  <- P           */
 46:   VecCopy(AP,ART);                 /*   ART <- AP          */
 47:   VecDotBegin(RT,ART,&btop);          /*   (RT,ART)           */

 49:   if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 50:     VecNormBegin(RT,NORM_2,&dp);        /*   dp <- RT'*RT       */
 51:     VecDotEnd   (RT,ART,&btop);           /*   (RT,ART)           */
 52:     VecNormEnd  (RT,NORM_2,&dp);        /*   dp <- RT'*RT       */
 53:   } else if (ksp->normtype == KSP_NORM_NONE) {
 54:       dp   = 0.0; /* meaningless value that is passed to monitor and convergence test */
 55:   } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 56:     VecNormBegin(R,NORM_2,&dp);         /*   dp <- R'*R         */
 57:     VecDotEnd   (RT,ART,&btop);          /*   (RT,ART)           */
 58:     VecNormEnd  (R,NORM_2,&dp);        /*   dp <- RT'*RT       */
 59:   } else if (ksp->normtype == KSP_NORM_NATURAL) {
 60:     VecDotEnd   (RT,ART,&btop);           /*   (RT,ART)           */
 61:     dp   = PetscSqrtReal(PetscAbsScalar(btop));                  /* dp = sqrt(R,AR)      */
 62:   } else SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype);
 63:   if (PetscAbsScalar(btop) < 0.0) {
 64:     ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
 65:     PetscInfo(ksp,"diverging due to indefinite or negative definite matrix\n");
 66:     return(0);
 67:   }

 69:   ksp->its   = 0;
 70:   KSPMonitor(ksp,0,dp);
 71:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
 72:   ksp->rnorm = dp;
 73:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
 74:   KSPLogResidualHistory(ksp,dp);
 75:   (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
 76:   if (ksp->reason) return(0);

 78:   i = 0;
 79:   do {
 80:     KSP_PCApply(ksp,AP,Q);  /*   Q <- B* AP          */

 82:     VecDot(AP,Q,&apq);
 83:     if (PetscRealPart(apq) <= 0.0) {
 84:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
 85:       PetscInfo(ksp,"KSPSolve_CR:diverging due to indefinite or negative definite PC\n");
 86:       break;
 87:     }
 88:     ai = btop/apq;                                      /* ai = (RT,ART)/(AP,Q)  */

 90:     VecAXPY(X,ai,P);              /*   X   <- X + ai*P     */
 91:     VecAXPY(RT,-ai,Q);             /*   RT  <- RT - ai*Q    */
 92:     KSP_MatMult(ksp,Amat,RT,ART);  /*   ART <-   A*RT       */
 93:     bbot = btop;
 94:     VecDotBegin(RT,ART,&btop);

 96:     if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 97:       VecNormBegin(RT,NORM_2,&dp);      /*   dp <- || RT ||      */
 98:       VecDotEnd   (RT,ART,&btop);
 99:       VecNormEnd  (RT,NORM_2,&dp);      /*   dp <- || RT ||      */
100:     } else if (ksp->normtype == KSP_NORM_NATURAL) {
101:       VecDotEnd(RT,ART,&btop);
102:       dp   = PetscSqrtReal(PetscAbsScalar(btop));                /* dp = sqrt(R,AR)       */
103:     } else if (ksp->normtype == KSP_NORM_NONE) {
104:       VecDotEnd(RT,ART,&btop);
105:       dp   = 0.0; /* meaningless value that is passed to monitor and convergence test */
106:     } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
107:       VecAXPY(R,ai,AP);           /*   R   <- R - ai*AP    */
108:       VecNormBegin(R,NORM_2,&dp);       /*   dp <- R'*R          */
109:       VecDotEnd   (RT,ART,&btop);
110:       VecNormEnd  (R,NORM_2,&dp);       /*   dp <- R'*R          */
111:     } else SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype);
112:     if (PetscAbsScalar(btop) < 0.0) {
113:       ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
114:       PetscInfo(ksp,"diverging due to indefinite or negative definite PC\n");
115:       break;
116:     }

118:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
119:     ksp->its++;
120:     ksp->rnorm = dp;
121:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

123:     KSPLogResidualHistory(ksp,dp);
124:     KSPMonitor(ksp,i+1,dp);
125:     (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
126:     if (ksp->reason) break;

128:     bi   = btop/bbot;
129:     VecAYPX(P,bi,RT);              /*   P <- RT + Bi P     */
130:     VecAYPX(AP,bi,ART);            /*   AP <- ART + Bi AP  */
131:     i++;
132:   } while (i<ksp->max_it);
133:   if (i >= ksp->max_it) ksp->reason =  KSP_DIVERGED_ITS;
134:   return(0);
135: }


138: /*MC
139:      KSPCR - This code implements the (preconditioned) conjugate residuals method

141:    Options Database Keys:
142: .   see KSPSolve()

144:    Level: beginner

146:    Notes: The operator and the preconditioner must be symmetric for this method. The
147:           preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE.
148:           Support only for left preconditioning.

150:    References:
151: .   1. - Magnus R. Hestenes and Eduard Stiefel, Methods of Conjugate Gradients for Solving Linear Systems, 
152:    Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379

154: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG
155: M*/
156: PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp)
157: {

161:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
162:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
163:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
164:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

166:   ksp->ops->setup          = KSPSetUp_CR;
167:   ksp->ops->solve          = KSPSolve_CR;
168:   ksp->ops->destroy        = KSPDestroyDefault;
169:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
170:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
171:   ksp->ops->setfromoptions = 0;
172:   ksp->ops->view           = 0;
173:   return(0);
174: }