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

  4: static PetscErrorCode KSPSetUp_BiCG(KSP ksp)
  5: {
  6:   KSPSetWorkVecs(ksp, 6);
  7:   return 0;
  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:   PCGetDiagonalScale(ksp->pc, &diagonalscale);

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

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

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

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

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

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

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

117:    Level: beginner

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

123:    Supports only left preconditioning

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

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

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

137:   ksp->ops->setup          = KSPSetUp_BiCG;
138:   ksp->ops->solve          = KSPSolve_BiCG;
139:   ksp->ops->destroy        = KSPDestroyDefault;
140:   ksp->ops->view           = NULL;
141:   ksp->ops->setfromoptions = NULL;
142:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
143:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
144:   return 0;
145: }