Actual source code: bicg.c

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

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

 10: static PetscErrorCode KSPSolve_BiCG(KSP ksp)
 11: {
 12:   PetscInt    i;
 13:   PetscBool   diagonalscale;
 14:   PetscScalar dpi, a = 1.0, beta, betaold = 1.0, b, ma;
 15:   PetscReal   dp;
 16:   Vec         X, B, Zl, Zr, Rl, Rr, Pl, Pr;
 17:   Mat         Amat, Pmat;

 19:   PetscFunctionBegin;
 20:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
 21:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

 23:   X  = ksp->vec_sol;
 24:   B  = ksp->vec_rhs;
 25:   Rl = ksp->work[0];
 26:   Zl = ksp->work[1];
 27:   Pl = ksp->work[2];
 28:   Rr = ksp->work[3];
 29:   Zr = ksp->work[4];
 30:   Pr = ksp->work[5];

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

 34:   if (!ksp->guess_zero) {
 35:     PetscCall(KSP_MatMult(ksp, Amat, X, Rr)); /*   r <- b - Ax       */
 36:     PetscCall(VecAYPX(Rr, -1.0, B));
 37:   } else {
 38:     PetscCall(VecCopy(B, Rr)); /*     r <- b (x is 0) */
 39:   }
 40:   PetscCall(VecCopy(Rr, Rl));
 41:   PetscCall(KSP_PCApply(ksp, Rr, Zr)); /*     z <- Br         */
 42:   PetscCall(KSP_PCApplyHermitianTranspose(ksp, Rl, Zl));
 43:   if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 44:     PetscCall(VecNorm(Zr, NORM_2, &dp)); /*    dp <- z'*z       */
 45:   } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 46:     PetscCall(VecNorm(Rr, NORM_2, &dp)); /*    dp <- r'*r       */
 47:   } else dp = 0.0;

 49:   KSPCheckNorm(ksp, dp);
 50:   PetscCall(KSPMonitor(ksp, 0, dp));
 51:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
 52:   ksp->its   = 0;
 53:   ksp->rnorm = dp;
 54:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
 55:   PetscCall(KSPLogResidualHistory(ksp, dp));
 56:   PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
 57:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

 59:   i = 0;
 60:   do {
 61:     PetscCall(VecDot(Zr, Rl, &beta)); /*     beta <- r'z     */
 62:     KSPCheckDot(ksp, beta);
 63:     if (!i) {
 64:       if (beta == 0.0) {
 65:         ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
 66:         PetscFunctionReturn(PETSC_SUCCESS);
 67:       }
 68:       PetscCall(VecCopy(Zr, Pr)); /*     p <- z          */
 69:       PetscCall(VecCopy(Zl, Pl));
 70:     } else {
 71:       b = beta / betaold;
 72:       PetscCall(VecAYPX(Pr, b, Zr)); /*     p <- z + b* p   */
 73:       b = PetscConj(b);
 74:       PetscCall(VecAYPX(Pl, b, Zl));
 75:     }
 76:     betaold = beta;
 77:     PetscCall(KSP_MatMult(ksp, Amat, Pr, Zr)); /*     z <- Kp         */
 78:     PetscCall(KSP_MatMultHermitianTranspose(ksp, Amat, Pl, Zl));
 79:     PetscCall(VecDot(Zr, Pl, &dpi)); /*     dpi <- z'p      */
 80:     KSPCheckDot(ksp, dpi);
 81:     a = beta / dpi;               /*     a = beta/p'z    */
 82:     PetscCall(VecAXPY(X, a, Pr)); /*     x <- x + ap     */
 83:     ma = -a;
 84:     PetscCall(VecAXPY(Rr, ma, Zr));
 85:     ma = PetscConj(ma);
 86:     PetscCall(VecAXPY(Rl, ma, Zl));
 87:     if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
 88:       PetscCall(KSP_PCApply(ksp, Rr, Zr)); /*     z <- Br         */
 89:       PetscCall(KSP_PCApplyHermitianTranspose(ksp, Rl, Zl));
 90:       PetscCall(VecNorm(Zr, NORM_2, &dp)); /*    dp <- z'*z       */
 91:     } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
 92:       PetscCall(VecNorm(Rr, NORM_2, &dp)); /*    dp <- r'*r       */
 93:     } else dp = 0.0;

 95:     KSPCheckNorm(ksp, dp);
 96:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
 97:     ksp->its   = i + 1;
 98:     ksp->rnorm = dp;
 99:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
100:     PetscCall(KSPLogResidualHistory(ksp, dp));
101:     PetscCall(KSPMonitor(ksp, i + 1, dp));
102:     PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
103:     if (ksp->reason) break;
104:     if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
105:       PetscCall(KSP_PCApply(ksp, Rr, Zr)); /* z <- Br  */
106:       PetscCall(KSP_PCApplyHermitianTranspose(ksp, Rl, Zl));
107:     }
108:     i++;
109:   } while (i < ksp->max_it);
110:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: /*MC
115:      KSPBICG - Implements the Biconjugate gradient method (similar to running the conjugate
116:          gradient on the normal equations).

118:    Level: beginner

120:    Notes:
121:    This method requires that one be apply to apply the transpose of the preconditioner and operator
122:    as well as the operator and preconditioner.

124:    Supports only left preconditioning

126:    See `KSPCGNE` for code that EXACTLY runs the preconditioned conjugate gradient method on the normal equations

128:    See `KSPBCGS` for the famous stabilized variant of this algorithm

130: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBCGS`, `KSPCGNE`
131: M*/
132: PETSC_EXTERN PetscErrorCode KSPCreate_BiCG(KSP ksp)
133: {
134:   PetscFunctionBegin;
135:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
136:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
137:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

139:   ksp->ops->setup          = KSPSetUp_BiCG;
140:   ksp->ops->solve          = KSPSolve_BiCG;
141:   ksp->ops->destroy        = KSPDestroyDefault;
142:   ksp->ops->view           = NULL;
143:   ksp->ops->setfromoptions = NULL;
144:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
145:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
146:   PetscFunctionReturn(PETSC_SUCCESS);
147: }