Actual source code: bicg.c


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

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

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

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

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

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

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

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

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

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

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

119:    Level: beginner

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

125:    Supports only left preconditioning

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

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

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

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