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