Actual source code: fbcgsr.c

      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>
     10: #include <petsc/private/vecimpl.h>

     12: static PetscErrorCode KSPSetUp_FBCGSR(KSP ksp)
     13: {

     17:   KSPSetWorkVecs(ksp,8);
     18:   return(0);
     19: }

     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;

     37:   if (!ksp->vec_rhs->petscnative) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Only coded for PETSc vectors");
     38:   VecGetLocalSize(ksp->vec_sol,&N);

     40:   X  = ksp->vec_sol;
     41:   B  = ksp->vec_rhs;
     42:   P2 = ksp->work[0];

     44:   /* The followings are involved in modified inner product calculations and vector updates */
     45:   RP = ksp->work[1]; VecGetArray(RP,(PetscScalar**)&rp); VecRestoreArray(RP,NULL);
     46:   R  = ksp->work[2]; VecGetArray(R,(PetscScalar**)&r);   VecRestoreArray(R,NULL);
     47:   P  = ksp->work[3]; VecGetArray(P,(PetscScalar**)&p);   VecRestoreArray(P,NULL);
     48:   V  = ksp->work[4]; VecGetArray(V,(PetscScalar**)&v);   VecRestoreArray(V,NULL);
     49:   S  = ksp->work[5]; VecGetArray(S,(PetscScalar**)&s);   VecRestoreArray(S,NULL);
     50:   T  = ksp->work[6]; VecGetArray(T,(PetscScalar**)&t);   VecRestoreArray(T,NULL);
     51:   S2 = ksp->work[7]; VecGetArray(S2,(PetscScalar**)&s2); VecRestoreArray(S2,NULL);

     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:       VecDuplicate(X,&bcgs->guess);
     58:     }
     59:     VecCopy(X,bcgs->guess);
     60:   } else {
     61:     VecSet(X,0.0);
     62:   }

     64:   /* Compute initial residual */
     65:   KSPGetPC(ksp,&pc);
     66:   PCSetUp(pc);
     67:   PCGetOperators(pc,&mat,NULL);
     68:   if (!ksp->guess_zero) {
     69:     KSP_MatMult(ksp,mat,X,P2); /* P2 is used as temporary storage */
     70:     VecCopy(B,R);
     71:     VecAXPY(R,-1.0,P2);
     72:   } else {
     73:     VecCopy(B,R);
     74:   }

     76:   /* Test for nothing to do */
     77:   VecNorm(R,NORM_2,&rho);
     78:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
     79:   ksp->its = 0;
     80:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rho;
     81:   else ksp->rnorm = 0;
     82:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
     83:   KSPLogResidualHistory(ksp,ksp->rnorm);
     84:   KSPMonitor(ksp,0,ksp->rnorm);
     85:   (*ksp->converged)(ksp,0,ksp->rnorm,&ksp->reason,ksp->cnvP);
     86:   if (ksp->reason) return(0);

     88:   /* Initialize iterates */
     89:   VecCopy(R,RP); /* rp <- r */
     90:   VecCopy(R,P); /* p <- r */

     92:   /* Big loop */
     93:   for (i=0; i<ksp->max_it; i++) {

     95:     /* matmult and pc */
     96:     KSP_PCApply(ksp,P,P2); /* p2 <- K p */
     97:     KSP_MatMult(ksp,mat,P2,V); /* v <- A p2 */

     99:     /* inner prodcuts */
    100:     if (i==0) {
    101:       tau  = rho*rho;
    102:       VecDot(V,RP,&sigma); /* sigma <- (v,rp) */
    103:     } else {
    104:       PetscLogEventBegin(VEC_ReduceArithmetic,0,0,0,0);
    105:       tau  = sigma = 0.0;
    106:       for (j=0; j<N; j++) {
    107:         tau   += r[j]*rp[j]; /* tau <- (r,rp) */
    108:         sigma += v[j]*rp[j]; /* sigma <- (v,rp) */
    109:       }
    110:       PetscLogFlops(4.0*N);
    111:       PetscLogEventEnd(VEC_ReduceArithmetic,0,0,0,0);
    112:       insums[0] = tau;
    113:       insums[1] = sigma;
    114:       PetscLogEventBegin(VEC_ReduceCommunication,0,0,0,0);
    115:       MPIU_Allreduce(insums,outsums,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
    116:       PetscLogEventEnd(VEC_ReduceCommunication,0,0,0,0);
    117:       tau       = outsums[0];
    118:       sigma     = outsums[1];
    119:     }

    121:     /* scalar update */
    122:     alpha = tau / sigma;

    124:     /* vector update */
    125:     VecWAXPY(S,-alpha,V,R);  /* s <- r - alpha v */

    127:     /* matmult and pc */
    128:     KSP_PCApply(ksp,S,S2); /* s2 <- K s */
    129:     KSP_MatMult(ksp,mat,S2,T); /* t <- A s2 */

    131:     /* inner prodcuts */
    132:     PetscLogEventBegin(VEC_ReduceArithmetic,0,0,0,0);
    133:     xi1  = xi2 = xi3 = xi4 = 0.0;
    134:     for (j=0; j<N; j++) {
    135:       xi1 += s[j]*s[j]; /* xi1 <- (s,s) */
    136:       xi2 += t[j]*s[j]; /* xi2 <- (t,s) */
    137:       xi3 += t[j]*t[j]; /* xi3 <- (t,t) */
    138:       xi4 += t[j]*rp[j]; /* xi4 <- (t,rp) */
    139:     }
    140:     PetscLogFlops(8.0*N);
    141:     PetscLogEventEnd(VEC_ReduceArithmetic,0,0,0,0);

    143:     insums[0] = xi1;
    144:     insums[1] = xi2;
    145:     insums[2] = xi3;
    146:     insums[3] = xi4;

    148:     PetscLogEventBegin(VEC_ReduceCommunication,0,0,0,0);
    149:     MPIU_Allreduce(insums,outsums,4,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
    150:     PetscLogEventEnd(VEC_ReduceCommunication,0,0,0,0);
    151:     xi1  = outsums[0];
    152:     xi2  = outsums[1];
    153:     xi3  = outsums[2];
    154:     xi4  = outsums[3];

    156:     /* test denominator */
    157:     if (xi3 == 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Divide by zero");
    158:     if (sigma == 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Divide by zero");

    160:     /* scalar updates */
    161:     omega = xi2 / xi3;
    162:     beta  = -xi4 / sigma;
    163:     rho   = PetscSqrtReal(PetscAbsScalar(xi1 - omega * xi2)); /* residual norm */

    165:     /* vector updates */
    166:     VecAXPBYPCZ(X,alpha,omega,1.0,P2,S2); /* x <- alpha * p2 + omega * s2 + x */

    168:     /* convergence test */
    169:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
    170:     ksp->its++;
    171:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rho;
    172:     else ksp->rnorm = 0;
    173:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
    174:     KSPLogResidualHistory(ksp,ksp->rnorm);
    175:     KSPMonitor(ksp,i+1,ksp->rnorm);
    176:     (*ksp->converged)(ksp,i+1,ksp->rnorm,&ksp->reason,ksp->cnvP);
    177:     if (ksp->reason) break;

    179:     /* vector updates */
    180:     PetscLogEventBegin(VEC_Ops,0,0,0,0);
    181:     for (j=0; j<N; j++) {
    182:       r[j] = s[j] - omega * t[j]; /* r <- s - omega t */
    183:       p[j] = r[j] + beta * (p[j] - omega * v[j]); /* p <- r + beta * (p - omega v) */
    184:     }
    185:     PetscLogFlops(6.0*N);
    186:     PetscLogEventEnd(VEC_Ops,0,0,0,0);

    188:   }

    190:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    191:   return(0);
    192: }

    194: /*MC
    195:      KSPFBCGSR - Implements a mathematically equivalent variant of FBiCGSTab.

    197:    Options Database Keys:
    198: .   see KSPSolve()

    200:    Level: beginner

    202:    Notes:
    203:     Only allow right preconditioning

    205: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPFBCGSL, KSPSetPCSide()
    206: M*/
    207: PETSC_EXTERN PetscErrorCode KSPCreate_FBCGSR(KSP ksp)
    208: {
    210:   KSP_BCGS       *bcgs;

    213:   PetscNewLog(ksp,&bcgs);

    215:   ksp->data                = bcgs;
    216:   ksp->ops->setup          = KSPSetUp_FBCGSR;
    217:   ksp->ops->solve          = KSPSolve_FBCGSR;
    218:   ksp->ops->destroy        = KSPDestroy_BCGS;
    219:   ksp->ops->reset          = KSPReset_BCGS;
    220:   ksp->ops->buildsolution  = KSPBuildSolution_BCGS;
    221:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    222:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
    223:   ksp->pc_side             = PC_RIGHT; /* set default PC side */

    225:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
    226:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
    227:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
    228:   return(0);
    229: }