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: }