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