Actual source code: bicg.c
2: #include <petsc/private/kspimpl.h>
4: static PetscErrorCode KSPSetUp_BiCG(KSP ksp)
5: {
6: /* check user parameters and functions */
9: KSPSetWorkVecs(ksp,6);
10: return 0;
11: }
13: static PetscErrorCode KSPSolve_BiCG(KSP ksp)
14: {
15: PetscInt i;
16: PetscBool diagonalscale;
17: PetscScalar dpi,a=1.0,beta,betaold=1.0,b,ma;
18: PetscReal dp;
19: Vec X,B,Zl,Zr,Rl,Rr,Pl,Pr;
20: Mat Amat,Pmat;
22: PCGetDiagonalScale(ksp->pc,&diagonalscale);
25: X = ksp->vec_sol;
26: B = ksp->vec_rhs;
27: Rl = ksp->work[0];
28: Zl = ksp->work[1];
29: Pl = ksp->work[2];
30: Rr = ksp->work[3];
31: Zr = ksp->work[4];
32: Pr = ksp->work[5];
34: PCGetOperators(ksp->pc,&Amat,&Pmat);
36: if (!ksp->guess_zero) {
37: KSP_MatMult(ksp,Amat,X,Rr); /* r <- b - Ax */
38: VecAYPX(Rr,-1.0,B);
39: } else {
40: VecCopy(B,Rr); /* r <- b (x is 0) */
41: }
42: VecCopy(Rr,Rl);
43: KSP_PCApply(ksp,Rr,Zr); /* z <- Br */
44: KSP_PCApplyHermitianTranspose(ksp,Rl,Zl);
45: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
46: VecNorm(Zr,NORM_2,&dp); /* dp <- z'*z */
47: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
48: VecNorm(Rr,NORM_2,&dp); /* dp <- r'*r */
49: } else dp = 0.0;
51: KSPCheckNorm(ksp,dp);
52: KSPMonitor(ksp,0,dp);
53: PetscObjectSAWsTakeAccess((PetscObject)ksp);
54: ksp->its = 0;
55: ksp->rnorm = dp;
56: PetscObjectSAWsGrantAccess((PetscObject)ksp);
57: KSPLogResidualHistory(ksp,dp);
58: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
59: if (ksp->reason) return 0;
61: i = 0;
62: do {
63: VecDot(Zr,Rl,&beta); /* beta <- r'z */
64: KSPCheckDot(ksp,beta);
65: if (!i) {
66: if (beta == 0.0) {
67: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
68: return 0;
69: }
70: VecCopy(Zr,Pr); /* p <- z */
71: VecCopy(Zl,Pl);
72: } else {
73: b = beta/betaold;
74: VecAYPX(Pr,b,Zr); /* p <- z + b* p */
75: b = PetscConj(b);
76: VecAYPX(Pl,b,Zl);
77: }
78: betaold = beta;
79: KSP_MatMult(ksp,Amat,Pr,Zr); /* z <- Kp */
80: KSP_MatMultHermitianTranspose(ksp,Amat,Pl,Zl);
81: VecDot(Zr,Pl,&dpi); /* dpi <- z'p */
82: KSPCheckDot(ksp,dpi);
83: a = beta/dpi; /* a = beta/p'z */
84: VecAXPY(X,a,Pr); /* x <- x + ap */
85: ma = -a;
86: VecAXPY(Rr,ma,Zr);
87: ma = PetscConj(ma);
88: VecAXPY(Rl,ma,Zl);
89: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
90: KSP_PCApply(ksp,Rr,Zr); /* z <- Br */
91: KSP_PCApplyHermitianTranspose(ksp,Rl,Zl);
92: VecNorm(Zr,NORM_2,&dp); /* dp <- z'*z */
93: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
94: VecNorm(Rr,NORM_2,&dp); /* dp <- r'*r */
95: } else dp = 0.0;
97: KSPCheckNorm(ksp,dp);
98: PetscObjectSAWsTakeAccess((PetscObject)ksp);
99: ksp->its = i+1;
100: ksp->rnorm = dp;
101: PetscObjectSAWsGrantAccess((PetscObject)ksp);
102: KSPLogResidualHistory(ksp,dp);
103: KSPMonitor(ksp,i+1,dp);
104: (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
105: if (ksp->reason) break;
106: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
107: KSP_PCApply(ksp,Rr,Zr); /* z <- Br */
108: KSP_PCApplyHermitianTranspose(ksp,Rl,Zl);
109: }
110: i++;
111: } while (i<ksp->max_it);
112: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
113: return 0;
114: }
116: /*MC
117: KSPBICG - Implements the Biconjugate gradient method (similar to running the conjugate
118: gradient on the normal equations).
120: Options Database Keys:
121: see KSPSolve()
123: Level: beginner
125: Notes:
126: this method requires that one be apply to apply the transpose of the preconditioner and operator
127: as well as the operator and preconditioner.
128: Supports only left preconditioning
130: See KSPCGNE for code that EXACTLY runs the preconditioned conjugate gradient method on the
131: normal equations
133: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBCGS, KSPCGNE
135: M*/
136: PETSC_EXTERN PetscErrorCode KSPCreate_BiCG(KSP ksp)
137: {
138: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
139: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
140: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
142: ksp->ops->setup = KSPSetUp_BiCG;
143: ksp->ops->solve = KSPSolve_BiCG;
144: ksp->ops->destroy = KSPDestroyDefault;
145: ksp->ops->view = NULL;
146: ksp->ops->setfromoptions = NULL;
147: ksp->ops->buildsolution = KSPBuildSolutionDefault;
148: ksp->ops->buildresidual = KSPBuildResidualDefault;
149: return 0;
150: }