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