Actual source code: bicg.c
petsc-3.9.4 2018-09-11
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: VecConjugate(Rl);
50: KSP_PCApplyTranspose(ksp,Rl,Zl);
51: VecConjugate(Rl);
52: VecConjugate(Zl);
53: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
54: VecNorm(Zr,NORM_2,&dp); /* dp <- z'*z */
55: } else {
56: VecNorm(Rr,NORM_2,&dp); /* dp <- r'*r */
57: }
58: KSPMonitor(ksp,0,dp);
59: PetscObjectSAWsTakeAccess((PetscObject)ksp);
60: ksp->its = 0;
61: ksp->rnorm = dp;
62: PetscObjectSAWsGrantAccess((PetscObject)ksp);
63: KSPLogResidualHistory(ksp,dp);
64: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
65: if (ksp->reason) return(0);
67: i = 0;
68: do {
69: VecDot(Zr,Rl,&beta); /* beta <- r'z */
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: VecConjugate(Pl);
86: KSP_MatMultTranspose(ksp,Amat,Pl,Zl);
87: VecConjugate(Pl);
88: VecConjugate(Zl);
89: VecDot(Zr,Pl,&dpi); /* dpi <- z'p */
90: a = beta/dpi; /* a = beta/p'z */
91: VecAXPY(X,a,Pr); /* x <- x + ap */
92: ma = -a;
93: VecAXPY(Rr,ma,Zr);
94: ma = PetscConj(ma);
95: VecAXPY(Rl,ma,Zl);
96: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
97: KSP_PCApply(ksp,Rr,Zr); /* z <- Br */
98: VecConjugate(Rl);
99: KSP_PCApplyTranspose(ksp,Rl,Zl);
100: VecConjugate(Rl);
101: VecConjugate(Zl);
102: VecNorm(Zr,NORM_2,&dp); /* dp <- z'*z */
103: } else {
104: VecNorm(Rr,NORM_2,&dp); /* dp <- r'*r */
105: }
106: PetscObjectSAWsTakeAccess((PetscObject)ksp);
107: ksp->its = i+1;
108: ksp->rnorm = dp;
109: PetscObjectSAWsGrantAccess((PetscObject)ksp);
110: KSPLogResidualHistory(ksp,dp);
111: KSPMonitor(ksp,i+1,dp);
112: (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
113: if (ksp->reason) break;
114: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
115: KSP_PCApply(ksp,Rr,Zr); /* z <- Br */
116: VecConjugate(Rl);
117: KSP_PCApplyTranspose(ksp,Rl,Zl);
118: VecConjugate(Rl);
119: VecConjugate(Zl);
120: }
121: i++;
122: } while (i<ksp->max_it);
123: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
124: return(0);
125: }
127: /*MC
128: KSPBICG - Implements the Biconjugate gradient method (similar to running the conjugate
129: gradient on the normal equations).
131: Options Database Keys:
132: . see KSPSolve()
134: Level: beginner
136: Notes: this method requires that one be apply to apply the transpose of the preconditioner and operator
137: as well as the operator and preconditioner.
138: Supports only left preconditioning
140: See KSPCGNE for code that EXACTLY runs the preconditioned conjugate gradient method on the
141: normal equations
143: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBCGS, KSPCGNE
145: M*/
146: PETSC_EXTERN PetscErrorCode KSPCreate_BiCG(KSP ksp)
147: {
151: ksp->data = (void*)0;
152: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
153: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
155: ksp->ops->setup = KSPSetUp_BiCG;
156: ksp->ops->solve = KSPSolve_BiCG;
157: ksp->ops->destroy = KSPDestroyDefault;
158: ksp->ops->view = 0;
159: ksp->ops->setfromoptions = 0;
160: ksp->ops->buildsolution = KSPBuildSolutionDefault;
161: ksp->ops->buildresidual = KSPBuildResidualDefault;
162: return(0);
163: }