1
2 #include <petsc/private/kspimpl.h>
3
KSPSetUp_BiCG(KSP ksp)4 static PetscErrorCode KSPSetUp_BiCG(KSP ksp)
5 {
6 PetscErrorCode ierr;
7
8 PetscFunctionBegin;
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 ierr = KSPSetWorkVecs(ksp,6);CHKERRQ(ierr);
13 PetscFunctionReturn(0);
14 }
15
KSPSolve_BiCG(KSP ksp)16 static PetscErrorCode KSPSolve_BiCG(KSP ksp)
17 {
18 PetscErrorCode ierr;
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;
25
26 PetscFunctionBegin;
27 ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
28 if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
29
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];
38
39 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
40
41 if (!ksp->guess_zero) {
42 ierr = KSP_MatMult(ksp,Amat,X,Rr);CHKERRQ(ierr); /* r <- b - Ax */
43 ierr = VecAYPX(Rr,-1.0,B);CHKERRQ(ierr);
44 } else {
45 ierr = VecCopy(B,Rr);CHKERRQ(ierr); /* r <- b (x is 0) */
46 }
47 ierr = VecCopy(Rr,Rl);CHKERRQ(ierr);
48 ierr = KSP_PCApply(ksp,Rr,Zr);CHKERRQ(ierr); /* z <- Br */
49 ierr = VecConjugate(Rl);CHKERRQ(ierr);
50 ierr = KSP_PCApplyTranspose(ksp,Rl,Zl);CHKERRQ(ierr);
51 ierr = VecConjugate(Rl);CHKERRQ(ierr);
52 ierr = VecConjugate(Zl);CHKERRQ(ierr);
53 if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
54 ierr = VecNorm(Zr,NORM_2,&dp);CHKERRQ(ierr); /* dp <- z'*z */
55 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
56 ierr = VecNorm(Rr,NORM_2,&dp);CHKERRQ(ierr); /* dp <- r'*r */
57 } else dp = 0.0;
58
59 KSPCheckNorm(ksp,dp);
60 ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr);
61 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
62 ksp->its = 0;
63 ksp->rnorm = dp;
64 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
65 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
66 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
67 if (ksp->reason) PetscFunctionReturn(0);
68
69 i = 0;
70 do {
71 ierr = VecDot(Zr,Rl,&beta);CHKERRQ(ierr); /* beta <- r'z */
72 KSPCheckDot(ksp,beta);
73 if (!i) {
74 if (beta == 0.0) {
75 ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
76 PetscFunctionReturn(0);
77 }
78 ierr = VecCopy(Zr,Pr);CHKERRQ(ierr); /* p <- z */
79 ierr = VecCopy(Zl,Pl);CHKERRQ(ierr);
80 } else {
81 b = beta/betaold;
82 ierr = VecAYPX(Pr,b,Zr);CHKERRQ(ierr); /* p <- z + b* p */
83 b = PetscConj(b);
84 ierr = VecAYPX(Pl,b,Zl);CHKERRQ(ierr);
85 }
86 betaold = beta;
87 ierr = KSP_MatMult(ksp,Amat,Pr,Zr);CHKERRQ(ierr); /* z <- Kp */
88 ierr = VecConjugate(Pl);CHKERRQ(ierr);
89 ierr = KSP_MatMultTranspose(ksp,Amat,Pl,Zl);CHKERRQ(ierr);
90 ierr = VecConjugate(Pl);CHKERRQ(ierr);
91 ierr = VecConjugate(Zl);CHKERRQ(ierr);
92 ierr = VecDot(Zr,Pl,&dpi);CHKERRQ(ierr); /* dpi <- z'p */
93 KSPCheckDot(ksp,dpi);
94 a = beta/dpi; /* a = beta/p'z */
95 ierr = VecAXPY(X,a,Pr);CHKERRQ(ierr); /* x <- x + ap */
96 ma = -a;
97 ierr = VecAXPY(Rr,ma,Zr);CHKERRQ(ierr);
98 ma = PetscConj(ma);
99 ierr = VecAXPY(Rl,ma,Zl);CHKERRQ(ierr);
100 if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
101 ierr = KSP_PCApply(ksp,Rr,Zr);CHKERRQ(ierr); /* z <- Br */
102 ierr = VecConjugate(Rl);CHKERRQ(ierr);
103 ierr = KSP_PCApplyTranspose(ksp,Rl,Zl);CHKERRQ(ierr);
104 ierr = VecConjugate(Rl);CHKERRQ(ierr);
105 ierr = VecConjugate(Zl);CHKERRQ(ierr);
106 ierr = VecNorm(Zr,NORM_2,&dp);CHKERRQ(ierr); /* dp <- z'*z */
107 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
108 ierr = VecNorm(Rr,NORM_2,&dp);CHKERRQ(ierr); /* dp <- r'*r */
109 } else dp = 0.0;
110
111 KSPCheckNorm(ksp,dp);
112 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
113 ksp->its = i+1;
114 ksp->rnorm = dp;
115 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
116 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr);
117 ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr);
118 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
119 if (ksp->reason) break;
120 if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
121 ierr = KSP_PCApply(ksp,Rr,Zr);CHKERRQ(ierr); /* z <- Br */
122 ierr = VecConjugate(Rl);CHKERRQ(ierr);
123 ierr = KSP_PCApplyTranspose(ksp,Rl,Zl);CHKERRQ(ierr);
124 ierr = VecConjugate(Rl);CHKERRQ(ierr);
125 ierr = VecConjugate(Zl);CHKERRQ(ierr);
126 }
127 i++;
128 } while (i<ksp->max_it);
129 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
130 PetscFunctionReturn(0);
131 }
132
133 /*MC
134 KSPBICG - Implements the Biconjugate gradient method (similar to running the conjugate
135 gradient on the normal equations).
136
137 Options Database Keys:
138 . see KSPSolve()
139
140 Level: beginner
141
142 Notes:
143 this method requires that one be apply to apply the transpose of the preconditioner and operator
144 as well as the operator and preconditioner.
145 Supports only left preconditioning
146
147 See KSPCGNE for code that EXACTLY runs the preconditioned conjugate gradient method on the
148 normal equations
149
150 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBCGS, KSPCGNE
151
152 M*/
KSPCreate_BiCG(KSP ksp)153 PETSC_EXTERN PetscErrorCode KSPCreate_BiCG(KSP ksp)
154 {
155 PetscErrorCode ierr;
156
157 PetscFunctionBegin;
158 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr);
159 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
160 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
161
162 ksp->ops->setup = KSPSetUp_BiCG;
163 ksp->ops->solve = KSPSolve_BiCG;
164 ksp->ops->destroy = KSPDestroyDefault;
165 ksp->ops->view = NULL;
166 ksp->ops->setfromoptions = NULL;
167 ksp->ops->buildsolution = KSPBuildSolutionDefault;
168 ksp->ops->buildresidual = KSPBuildResidualDefault;
169 PetscFunctionReturn(0);
170 }
171