Actual source code: cr.c


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

  4: static PetscErrorCode KSPSetUp_CR(KSP ksp)
  5: {
  6:   PetscFunctionBegin;
  7:   PetscCall(KSPSetWorkVecs(ksp, 6));
  8:   PetscFunctionReturn(PETSC_SUCCESS);
  9: }

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

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

 30:   /* R is the true residual norm, RT is the preconditioned residual norm */
 31:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
 32:   if (!ksp->guess_zero) {
 33:     PetscCall(KSP_MatMult(ksp, Amat, X, R)); /*   R <- A*X           */
 34:     PetscCall(VecAYPX(R, -1.0, B));          /*   R <- B-R == B-A*X  */
 35:   } else {
 36:     PetscCall(VecCopy(B, R)); /*   R <- B (X is 0)    */
 37:   }
 38:   /* 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 */
 39:   if (ksp->reason == KSP_DIVERGED_PC_FAILED) PetscCall(VecSetInf(R));
 40:   PetscCall(KSP_PCApply(ksp, R, P));        /*   P   <- B*R         */
 41:   PetscCall(KSP_MatMult(ksp, Amat, P, AP)); /*   AP  <- A*P         */
 42:   PetscCall(VecCopy(P, RT));                /*   RT  <- P           */
 43:   PetscCall(VecCopy(AP, ART));              /*   ART <- AP          */
 44:   PetscCall(VecDotBegin(RT, ART, &btop));   /*   (RT,ART)           */

 46:   if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 47:     PetscCall(VecNormBegin(RT, NORM_2, &dp)); /*   dp <- RT'*RT       */
 48:     PetscCall(VecDotEnd(RT, ART, &btop));     /*   (RT,ART)           */
 49:     PetscCall(VecNormEnd(RT, NORM_2, &dp));   /*   dp <- RT'*RT       */
 50:     KSPCheckNorm(ksp, dp);
 51:   } else if (ksp->normtype == KSP_NORM_NONE) {
 52:     dp = 0.0;                             /* meaningless value that is passed to monitor and convergence test */
 53:     PetscCall(VecDotEnd(RT, ART, &btop)); /*   (RT,ART)           */
 54:   } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 55:     PetscCall(VecNormBegin(R, NORM_2, &dp)); /*   dp <- R'*R         */
 56:     PetscCall(VecDotEnd(RT, ART, &btop));    /*   (RT,ART)           */
 57:     PetscCall(VecNormEnd(R, NORM_2, &dp));   /*   dp <- RT'*RT       */
 58:     KSPCheckNorm(ksp, dp);
 59:   } else if (ksp->normtype == KSP_NORM_NATURAL) {
 60:     PetscCall(VecDotEnd(RT, ART, &btop));     /*   (RT,ART)           */
 61:     dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR)      */
 62:   } else SETERRQ(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:     PetscCall(PetscInfo(ksp, "diverging due to indefinite or negative definite matrix\n"));
 66:     PetscFunctionReturn(PETSC_SUCCESS);
 67:   }

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

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

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

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

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

121:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
122:     ksp->its++;
123:     ksp->rnorm = dp;
124:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

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

131:     bi = btop / bbot;
132:     PetscCall(VecAYPX(P, bi, RT));   /*   P <- RT + Bi P     */
133:     PetscCall(VecAYPX(AP, bi, ART)); /*   AP <- ART + Bi AP  */
134:     i++;
135:   } while (i < ksp->max_it);
136:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

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

143:    Level: beginner

145:    Notes:
146:    The operator and the preconditioner must be symmetric for this method.

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

150:    Support only for left preconditioning.

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

156: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`
157: M*/
158: PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp)
159: {
160:   PetscFunctionBegin;
161:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
162:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
163:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
164:   PetscCall(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 = NULL;
172:   ksp->ops->view           = NULL;
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }