1 #include <petsc/private/kspimpl.h>
2
KSPSetUp_BiCG(KSP ksp)3 static PetscErrorCode KSPSetUp_BiCG(KSP ksp)
4 {
5 PetscFunctionBegin;
6 PetscCall(KSPSetWorkVecs(ksp, 6));
7 PetscFunctionReturn(PETSC_SUCCESS);
8 }
9
KSPSolve_BiCG(KSP ksp)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;
18
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);
22
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];
31
32 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
33
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;
48
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);
58
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;
94
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 }
113
114 /*MC
115 KSPBICG - Implements the Biconjugate gradient method (similar to running the conjugate
116 gradient on the normal equations).
117
118 Level: beginner
119
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.
123
124 Supports only left preconditioning
125
126 See `KSPCGNE` for code that EXACTLY runs the preconditioned conjugate gradient method on the normal equations
127
128 See `KSPBCGS` for the famous stabilized variant of this algorithm
129
130 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBCGS`, `KSPCGNE`
131 M*/
KSPCreate_BiCG(KSP ksp)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));
138
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 }
148