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