1 2 /* 3 This file implements FBiCGStab-R. 4 Only allow right preconditioning. 5 FBiCGStab-R is a mathematically equivalent variant of FBiCGStab. Differences are: 6 (1) There are fewer MPI_Allreduce calls. 7 (2) The convergence occasionally is much faster than that of FBiCGStab. 8 */ 9 #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h> /*I "petscksp.h" I*/ 10 #include <petsc/private/vecimpl.h> 11 12 static PetscErrorCode KSPSetUp_FBCGSR(KSP ksp) 13 { 14 PetscErrorCode ierr; 15 16 PetscFunctionBegin; 17 ierr = KSPSetWorkVecs(ksp,8);CHKERRQ(ierr); 18 PetscFunctionReturn(0); 19 } 20 21 static PetscErrorCode KSPSolve_FBCGSR(KSP ksp) 22 { 23 PetscErrorCode ierr; 24 PetscInt i,j,N; 25 PetscScalar tau,sigma,alpha,omega,beta; 26 PetscReal rho; 27 PetscScalar xi1,xi2,xi3,xi4; 28 Vec X,B,P,P2,RP,R,V,S,T,S2; 29 PetscScalar *PETSC_RESTRICT rp, *PETSC_RESTRICT r, *PETSC_RESTRICT p; 30 PetscScalar *PETSC_RESTRICT v, *PETSC_RESTRICT s, *PETSC_RESTRICT t, *PETSC_RESTRICT s2; 31 PetscScalar insums[4],outsums[4]; 32 KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data; 33 PC pc; 34 Mat mat; 35 36 PetscFunctionBegin; 37 if (!ksp->vec_rhs->petscnative) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Only coded for PETSc vectors"); 38 ierr = VecGetLocalSize(ksp->vec_sol,&N);CHKERRQ(ierr); 39 40 X = ksp->vec_sol; 41 B = ksp->vec_rhs; 42 P2 = ksp->work[0]; 43 44 /* The followings are involved in modified inner product calculations and vector updates */ 45 RP = ksp->work[1]; ierr = VecGetArray(RP,(PetscScalar**)&rp);CHKERRQ(ierr); ierr = VecRestoreArray(RP,NULL);CHKERRQ(ierr); 46 R = ksp->work[2]; ierr = VecGetArray(R,(PetscScalar**)&r);CHKERRQ(ierr); ierr = VecRestoreArray(R,NULL);CHKERRQ(ierr); 47 P = ksp->work[3]; ierr = VecGetArray(P,(PetscScalar**)&p);CHKERRQ(ierr); ierr = VecRestoreArray(P,NULL);CHKERRQ(ierr); 48 V = ksp->work[4]; ierr = VecGetArray(V,(PetscScalar**)&v);CHKERRQ(ierr); ierr = VecRestoreArray(V,NULL);CHKERRQ(ierr); 49 S = ksp->work[5]; ierr = VecGetArray(S,(PetscScalar**)&s);CHKERRQ(ierr); ierr = VecRestoreArray(S,NULL);CHKERRQ(ierr); 50 T = ksp->work[6]; ierr = VecGetArray(T,(PetscScalar**)&t);CHKERRQ(ierr); ierr = VecRestoreArray(T,NULL);CHKERRQ(ierr); 51 S2 = ksp->work[7]; ierr = VecGetArray(S2,(PetscScalar**)&s2);CHKERRQ(ierr); ierr = VecRestoreArray(S2,NULL);CHKERRQ(ierr); 52 53 /* Only supports right preconditioning */ 54 if (ksp->pc_side != PC_RIGHT) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP fbcgsr does not support %s",PCSides[ksp->pc_side]); 55 if (!ksp->guess_zero) { 56 if (!bcgs->guess) { 57 ierr = VecDuplicate(X,&bcgs->guess);CHKERRQ(ierr); 58 } 59 ierr = VecCopy(X,bcgs->guess);CHKERRQ(ierr); 60 } else { 61 ierr = VecSet(X,0.0);CHKERRQ(ierr); 62 } 63 64 /* Compute initial residual */ 65 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 66 ierr = PCSetUp(pc);CHKERRQ(ierr); 67 ierr = PCGetOperators(pc,&mat,NULL);CHKERRQ(ierr); 68 if (!ksp->guess_zero) { 69 ierr = KSP_MatMult(ksp,mat,X,P2);CHKERRQ(ierr); /* P2 is used as temporary storage */ 70 ierr = VecCopy(B,R);CHKERRQ(ierr); 71 ierr = VecAXPY(R,-1.0,P2);CHKERRQ(ierr); 72 } else { 73 ierr = VecCopy(B,R);CHKERRQ(ierr); 74 } 75 76 /* Test for nothing to do */ 77 ierr = VecNorm(R,NORM_2,&rho);CHKERRQ(ierr); 78 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 79 ksp->its = 0; 80 ksp->rnorm = rho; 81 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 82 ierr = KSPLogResidualHistory(ksp,rho);CHKERRQ(ierr); 83 ierr = KSPMonitor(ksp,0,rho);CHKERRQ(ierr); 84 ierr = (*ksp->converged)(ksp,0,rho,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 85 if (ksp->reason) PetscFunctionReturn(0); 86 87 /* Initialize iterates */ 88 ierr = VecCopy(R,RP);CHKERRQ(ierr); /* rp <- r */ 89 ierr = VecCopy(R,P);CHKERRQ(ierr); /* p <- r */ 90 91 /* Big loop */ 92 for (i=0; i<ksp->max_it; i++) { 93 94 /* matmult and pc */ 95 ierr = KSP_PCApply(ksp,P,P2);CHKERRQ(ierr); /* p2 <- K p */ 96 ierr = KSP_MatMult(ksp,mat,P2,V);CHKERRQ(ierr); /* v <- A p2 */ 97 98 /* inner prodcuts */ 99 if (i==0) { 100 tau = rho*rho; 101 ierr = VecDot(V,RP,&sigma);CHKERRQ(ierr); /* sigma <- (v,rp) */ 102 } else { 103 ierr = PetscLogEventBegin(VEC_ReduceArithmetic,0,0,0,0);CHKERRQ(ierr); 104 tau = sigma = 0.0; 105 for (j=0; j<N; j++) { 106 tau += r[j]*rp[j]; /* tau <- (r,rp) */ 107 sigma += v[j]*rp[j]; /* sigma <- (v,rp) */ 108 } 109 ierr = PetscLogFlops(4.0*N);CHKERRQ(ierr); 110 ierr = PetscLogEventEnd(VEC_ReduceArithmetic,0,0,0,0);CHKERRQ(ierr); 111 insums[0] = tau; 112 insums[1] = sigma; 113 ierr = PetscLogEventBarrierBegin(VEC_ReduceBarrier,0,0,0,0,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); 114 ierr = MPIU_Allreduce(insums,outsums,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); 115 ierr = PetscLogEventBarrierEnd(VEC_ReduceBarrier,0,0,0,0,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); 116 tau = outsums[0]; 117 sigma = outsums[1]; 118 } 119 120 /* scalar update */ 121 alpha = tau / sigma; 122 123 /* vector update */ 124 ierr = VecWAXPY(S,-alpha,V,R);CHKERRQ(ierr); /* s <- r - alpha v */ 125 126 /* matmult and pc */ 127 ierr = KSP_PCApply(ksp,S,S2);CHKERRQ(ierr); /* s2 <- K s */ 128 ierr = KSP_MatMult(ksp,mat,S2,T);CHKERRQ(ierr); /* t <- A s2 */ 129 130 /* inner prodcuts */ 131 ierr = PetscLogEventBegin(VEC_ReduceArithmetic,0,0,0,0);CHKERRQ(ierr); 132 xi1 = xi2 = xi3 = xi4 = 0.0; 133 for (j=0; j<N; j++) { 134 xi1 += s[j]*s[j]; /* xi1 <- (s,s) */ 135 xi2 += t[j]*s[j]; /* xi2 <- (t,s) */ 136 xi3 += t[j]*t[j]; /* xi3 <- (t,t) */ 137 xi4 += t[j]*rp[j]; /* xi4 <- (t,rp) */ 138 } 139 ierr = PetscLogFlops(8.0*N);CHKERRQ(ierr); 140 ierr = PetscLogEventEnd(VEC_ReduceArithmetic,0,0,0,0);CHKERRQ(ierr); 141 142 insums[0] = xi1; 143 insums[1] = xi2; 144 insums[2] = xi3; 145 insums[3] = xi4; 146 147 ierr = PetscLogEventBarrierBegin(VEC_ReduceBarrier,0,0,0,0,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); 148 ierr = MPIU_Allreduce(insums,outsums,4,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); 149 ierr = PetscLogEventBarrierEnd(VEC_ReduceBarrier,0,0,0,0,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); 150 xi1 = outsums[0]; 151 xi2 = outsums[1]; 152 xi3 = outsums[2]; 153 xi4 = outsums[3]; 154 155 /* test denominator */ 156 if (xi3 == 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Divide by zero"); 157 if (sigma == 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Divide by zero"); 158 159 /* scalar updates */ 160 omega = xi2 / xi3; 161 beta = -xi4 / sigma; 162 rho = PetscSqrtReal(PetscAbsScalar(xi1 - omega * xi2)); /* residual norm */ 163 164 /* vector updates */ 165 ierr = VecAXPBYPCZ(X,alpha,omega,1.0,P2,S2);CHKERRQ(ierr); /* x <- alpha * p2 + omega * s2 + x */ 166 167 /* convergence test */ 168 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 169 ksp->its++; 170 ksp->rnorm = rho; 171 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 172 ierr = KSPLogResidualHistory(ksp,rho);CHKERRQ(ierr); 173 ierr = KSPMonitor(ksp,i+1,rho);CHKERRQ(ierr); 174 ierr = (*ksp->converged)(ksp,i+1,rho,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 175 if (ksp->reason) break; 176 177 /* vector updates */ 178 ierr = PetscLogEventBegin(VEC_Ops,0,0,0,0);CHKERRQ(ierr); 179 for (j=0; j<N; j++) { 180 r[j] = s[j] - omega * t[j]; /* r <- s - omega t */ 181 p[j] = r[j] + beta * (p[j] - omega * v[j]); /* p <- r + beta * (p - omega v) */ 182 } 183 ierr = PetscLogFlops(6.0*N);CHKERRQ(ierr); 184 ierr = PetscLogEventEnd(VEC_Ops,0,0,0,0);CHKERRQ(ierr); 185 186 } 187 188 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 189 PetscFunctionReturn(0); 190 } 191 192 /*MC 193 KSPFBCGSR - Implements a mathematically equivalent variant of FBiCGSTab. 194 195 Options Database Keys: 196 . see KSPSolve() 197 198 Level: beginner 199 200 Notes: Only allow right preconditioning 201 202 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPFBCGSL, KSPSetPCSide() 203 M*/ 204 PETSC_EXTERN PetscErrorCode KSPCreate_FBCGSR(KSP ksp) 205 { 206 PetscErrorCode ierr; 207 KSP_BCGS *bcgs; 208 209 PetscFunctionBegin; 210 ierr = PetscNewLog(ksp,&bcgs);CHKERRQ(ierr); 211 212 ksp->data = bcgs; 213 ksp->ops->setup = KSPSetUp_FBCGSR; 214 ksp->ops->solve = KSPSolve_FBCGSR; 215 ksp->ops->destroy = KSPDestroy_BCGS; 216 ksp->ops->reset = KSPReset_BCGS; 217 ksp->ops->buildsolution = KSPBuildSolution_BCGS; 218 ksp->ops->buildresidual = KSPBuildResidualDefault; 219 ksp->ops->setfromoptions = KSPSetFromOptions_BCGS; 220 ksp->pc_side = PC_RIGHT; /* set default PC side */ 221 222 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr); 223 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr); 224 PetscFunctionReturn(0); 225 } 226