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