1 2 /* 3 This file implements flexible BiCGStab (FBiCGStab). 4 Only allow right preconditioning. 5 */ 6 #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h> /*I "petscksp.h" I*/ 7 8 /* copied from KSPBuildSolution_GCR() */ 9 #undef __FUNCT__ 10 #define __FUNCT__ "KSPBuildSolution_FBCGS" 11 PetscErrorCode KSPBuildSolution_FBCGS(KSP ksp, Vec v, Vec *V) 12 { 13 PetscErrorCode ierr; 14 Vec x; 15 16 PetscFunctionBegin; 17 x = ksp->vec_sol; 18 if (v) { 19 ierr = VecCopy(x, v);CHKERRQ(ierr); 20 if (V) *V = v; 21 } else if (V) { 22 *V = ksp->vec_sol; 23 } 24 PetscFunctionReturn(0); 25 } 26 27 #undef __FUNCT__ 28 #define __FUNCT__ "KSPSetUp_FBCGS" 29 static PetscErrorCode KSPSetUp_FBCGS(KSP ksp) 30 { 31 PetscErrorCode ierr; 32 33 PetscFunctionBegin; 34 ierr = KSPDefaultGetWork(ksp,8);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 /* Only need a few hacks from KSPSolve_BCGS */ 39 #include <petsc-private/pcimpl.h> /*I "petscksp.h" I*/ 40 #undef __FUNCT__ 41 #define __FUNCT__ "KSPSolve_FBCGS" 42 static PetscErrorCode KSPSolve_FBCGS(KSP ksp) 43 { 44 PetscErrorCode ierr; 45 PetscInt i; 46 PetscScalar rho,rhoold,alpha,beta,omega,omegaold,d1; 47 Vec X,B,V,P,R,RP,T,S,P2,S2; 48 PetscReal dp = 0.0,d2; 49 KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data; 50 PC pc; 51 52 PetscFunctionBegin; 53 X = ksp->vec_sol; 54 B = ksp->vec_rhs; 55 R = ksp->work[0]; 56 RP = ksp->work[1]; 57 V = ksp->work[2]; 58 T = ksp->work[3]; 59 S = ksp->work[4]; 60 P = ksp->work[5]; 61 S2 = ksp->work[6]; 62 P2 = ksp->work[7]; 63 64 /* Only supports right preconditioning */ 65 if (ksp->pc_side != PC_RIGHT) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"KSP fbcgs does not support %s",PCSides[ksp->pc_side]); 66 if (!ksp->guess_zero) { 67 if (!bcgs->guess) { 68 ierr = VecDuplicate(X,&bcgs->guess);CHKERRQ(ierr); 69 } 70 ierr = VecCopy(X,bcgs->guess);CHKERRQ(ierr); 71 } else { 72 ierr = VecSet(X,0.0);CHKERRQ(ierr); 73 } 74 75 /* Compute initial residual */ 76 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 77 ierr = PCSetUp(pc);CHKERRQ(ierr); 78 if (!ksp->guess_zero) { 79 ierr = MatMult(pc->mat,X,S2);CHKERRQ(ierr); 80 ierr = VecCopy(B,R);CHKERRQ(ierr); 81 ierr = VecAXPY(R,-1.0,S2);CHKERRQ(ierr); 82 } else { 83 ierr = VecCopy(B,R);CHKERRQ(ierr); 84 } 85 86 /* Test for nothing to do */ 87 if (ksp->normtype != KSP_NORM_NONE) { 88 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); 89 } 90 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 91 ksp->its = 0; 92 ksp->rnorm = dp; 93 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 94 KSPLogResidualHistory(ksp,dp); 95 ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr); 96 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 97 if (ksp->reason) PetscFunctionReturn(0); 98 99 /* Make the initial Rp == R */ 100 ierr = VecCopy(R,RP);CHKERRQ(ierr); 101 102 rhoold = 1.0; 103 alpha = 1.0; 104 omegaold = 1.0; 105 ierr = VecSet(P,0.0);CHKERRQ(ierr); 106 ierr = VecSet(V,0.0);CHKERRQ(ierr); 107 108 i=0; 109 do { 110 ierr = VecDot(R,RP,&rho);CHKERRQ(ierr); /* rho <- (r,rp) */ 111 beta = (rho/rhoold) * (alpha/omegaold); 112 ierr = VecAXPBYPCZ(P,1.0,-omegaold*beta,beta,R,V);CHKERRQ(ierr); /* p <- r - omega * beta* v + beta * p */ 113 114 ierr = PCApply(pc,P,P2);CHKERRQ(ierr); /* p2 <- K p */ 115 ierr = MatMult(pc->mat,P2,V);CHKERRQ(ierr); /* v <- A p2 */ 116 117 ierr = VecDot(V,RP,&d1);CHKERRQ(ierr); 118 if (d1 == 0.0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_PLIB,"Divide by zero"); 119 alpha = rho / d1; /* alpha <- rho / (v,rp) */ 120 ierr = VecWAXPY(S,-alpha,V,R);CHKERRQ(ierr); /* s <- r - alpha v */ 121 122 ierr = PCApply(pc,S,S2);CHKERRQ(ierr); /* s2 <- K s */ 123 ierr = MatMult(pc->mat,S2,T);CHKERRQ(ierr); /* t <- A s2 */ 124 125 ierr = VecDotNorm2(S,T,&d1,&d2);CHKERRQ(ierr); 126 if (d2 == 0.0) { 127 /* t is 0. if s is 0, then alpha v == r, and hence alpha p may be our solution. Give it a try? */ 128 ierr = VecDot(S,S,&d1);CHKERRQ(ierr); 129 if (d1 != 0.0) { 130 ksp->reason = KSP_DIVERGED_BREAKDOWN; 131 break; 132 } 133 ierr = VecAXPY(X,alpha,P2);CHKERRQ(ierr); /* x <- x + alpha p2 */ 134 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 135 ksp->its++; 136 ksp->rnorm = 0.0; 137 ksp->reason = KSP_CONVERGED_RTOL; 138 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 139 KSPLogResidualHistory(ksp,dp); 140 ierr = KSPMonitor(ksp,i+1,0.0);CHKERRQ(ierr); 141 break; 142 } 143 omega = d1 / d2; /* omega <- (t's) / (t't) */ 144 ierr = VecAXPBYPCZ(X,alpha,omega,1.0,P2,S2);CHKERRQ(ierr); /* x <- alpha * p2 + omega * s2 + x */ 145 146 ierr = VecWAXPY(R,-omega,T,S);CHKERRQ(ierr); /* r <- s - omega t */ 147 if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i+2) { 148 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); 149 } 150 151 rhoold = rho; 152 omegaold = omega; 153 154 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 155 ksp->its++; 156 ksp->rnorm = dp; 157 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 158 KSPLogResidualHistory(ksp,dp); 159 ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr); 160 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 161 if (ksp->reason) break; 162 if (rho == 0.0) { 163 ksp->reason = KSP_DIVERGED_BREAKDOWN; 164 break; 165 } 166 i++; 167 } while (i<ksp->max_it); 168 169 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 170 PetscFunctionReturn(0); 171 } 172 173 /*MC 174 KSPFBCGS - Implements flexible BiCGStab method. 175 176 Options Database Keys: 177 . see KSPSolve() 178 179 Level: beginner 180 181 Notes: Only allow right preconditioning 182 183 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPFBCGSL, KSPSetPCSide() 184 M*/ 185 EXTERN_C_BEGIN 186 #undef __FUNCT__ 187 #define __FUNCT__ "KSPCreate_FBCGS" 188 PetscErrorCode KSPCreate_FBCGS(KSP ksp) 189 { 190 PetscErrorCode ierr; 191 KSP_BCGS *bcgs; 192 193 PetscFunctionBegin; 194 ierr = PetscNewLog(ksp,KSP_BCGS,&bcgs);CHKERRQ(ierr); 195 196 ksp->data = bcgs; 197 ksp->ops->setup = KSPSetUp_FBCGS; 198 ksp->ops->solve = KSPSolve_FBCGS; 199 ksp->ops->destroy = KSPDestroy_BCGS; 200 ksp->ops->reset = KSPReset_BCGS; 201 ksp->ops->buildsolution = KSPBuildSolution_FBCGS; 202 ksp->ops->buildresidual = KSPDefaultBuildResidual; 203 ksp->ops->setfromoptions = KSPSetFromOptions_BCGS; 204 ksp->pc_side = PC_RIGHT; /* set default PC side */ 205 206 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 207 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);CHKERRQ(ierr); 208 PetscFunctionReturn(0); 209 } 210 EXTERN_C_END 211