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:   PetscScalar dpi, a = 1.0, beta, betaold = 1.0, b, ma;
 14:   PetscReal   dp;
 15:   Vec         X, B, Zl, Zr, Rl, Rr, Pl, Pr;
 16:   Mat         Amat, Pmat;

 18:   PetscFunctionBegin;
 19:   X  = ksp->vec_sol;
 20:   B  = ksp->vec_rhs;
 21:   Rl = ksp->work[0];
 22:   Zl = ksp->work[1];
 23:   Pl = ksp->work[2];
 24:   Rr = ksp->work[3];
 25:   Zr = ksp->work[4];
 26:   Pr = ksp->work[5];

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

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

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

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

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

110: /*MC
111:      KSPBICG - Implements the Biconjugate gradient method (similar to running the conjugate
112:          gradient on the normal equations).

114:    Level: beginner

116:    Notes:
117:    This method requires that one be apply to apply the transpose of the preconditioner and operator
118:    as well as the operator and preconditioner.

120:    Supports only left preconditioning

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

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

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

135:   ksp->ops->setup          = KSPSetUp_BiCG;
136:   ksp->ops->solve          = KSPSolve_BiCG;
137:   ksp->ops->destroy        = KSPDestroyDefault;
138:   ksp->ops->view           = NULL;
139:   ksp->ops->setfromoptions = NULL;
140:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
141:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }