Actual source code: cr.c


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

  4: static PetscErrorCode KSPSetUp_CR(KSP ksp)
  5: {
  6:   KSPSetWorkVecs(ksp, 6);
  7:   return 0;
  8: }

 10: static PetscErrorCode KSPSolve_CR(KSP ksp)
 11: {
 12:   PetscInt    i = 0;
 13:   PetscReal   dp;
 14:   PetscScalar ai, bi;
 15:   PetscScalar apq, btop, bbot;
 16:   Vec         X, B, R, RT, P, AP, ART, Q;
 17:   Mat         Amat, Pmat;

 19:   X   = ksp->vec_sol;
 20:   B   = ksp->vec_rhs;
 21:   R   = ksp->work[0];
 22:   RT  = ksp->work[1];
 23:   P   = ksp->work[2];
 24:   AP  = ksp->work[3];
 25:   ART = ksp->work[4];
 26:   Q   = ksp->work[5];

 28:   /* R is the true residual norm, RT is the preconditioned residual norm */
 29:   PCGetOperators(ksp->pc, &Amat, &Pmat);
 30:   if (!ksp->guess_zero) {
 31:     KSP_MatMult(ksp, Amat, X, R); /*   R <- A*X           */
 32:     VecAYPX(R, -1.0, B);          /*   R <- B-R == B-A*X  */
 33:   } else {
 34:     VecCopy(B, R); /*   R <- B (X is 0)    */
 35:   }
 36:   /* This may be true only on a subset of MPI ranks; setting it here so it will be detected by the first norm computation below */
 37:   if (ksp->reason == KSP_DIVERGED_PC_FAILED) VecSetInf(R);
 38:   KSP_PCApply(ksp, R, P);        /*   P   <- B*R         */
 39:   KSP_MatMult(ksp, Amat, P, AP); /*   AP  <- A*P         */
 40:   VecCopy(P, RT);                /*   RT  <- P           */
 41:   VecCopy(AP, ART);              /*   ART <- AP          */
 42:   VecDotBegin(RT, ART, &btop);   /*   (RT,ART)           */

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

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

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

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

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

 95:     if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 96:       VecNormBegin(RT, NORM_2, &dp); /*   dp <- || RT ||      */
 97:       VecDotEnd(RT, ART, &btop);
 98:       VecNormEnd(RT, NORM_2, &dp); /*   dp <- || RT ||      */
 99:       KSPCheckNorm(ksp, dp);
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:       KSPCheckNorm(ksp, dp);
112:     } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSPNormType of %d not supported", (int)ksp->normtype);
113:     if (PetscAbsScalar(btop) < 0.0) {
114:       ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
115:       PetscInfo(ksp, "diverging due to indefinite or negative definite PC\n");
116:       break;
117:     }

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

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

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

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

141:    Level: beginner

143:    Notes:
144:    The operator and the preconditioner must be symmetric for this method.

146:    The preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE.

148:    Support only for left preconditioning.

150:    References:
151: .  * - 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: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`
155: M*/
156: PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp)
157: {
158:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
159:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);
160:   KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
161:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);

163:   ksp->ops->setup          = KSPSetUp_CR;
164:   ksp->ops->solve          = KSPSolve_CR;
165:   ksp->ops->destroy        = KSPDestroyDefault;
166:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
167:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
168:   ksp->ops->setfromoptions = NULL;
169:   ksp->ops->view           = NULL;
170:   return 0;
171: }