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