Actual source code: lsqr.c

      1: #define PETSCKSP_DLL

      3: #define SWAP(a,b,c) { c = a; a = b; b = c; }

     5:  #include src/ksp/ksp/kspimpl.h

      7: typedef struct {
      8:   PetscInt  nwork_n,nwork_m;
      9:   Vec       *vwork_m;  /* work vectors of length m, where the system is size m x n */
     10:   Vec       *vwork_n;  /* work vectors of length m */
     11: } KSP_LSQR;

     15: static PetscErrorCode KSPSetUp_LSQR(KSP ksp)
     16: {
     18:   PetscInt       nw;
     19:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

     22:   if (ksp->pc_side == PC_SYMMETRIC){
     23:     SETERRQ(PETSC_ERR_SUP,"no symmetric preconditioning for KSPLSQR");
     24:   }

     26:   /* Get work vectors */
     27:   lsqr->nwork_m = nw = 2;
     28:   if (lsqr->vwork_m) {
     29:     VecDestroyVecs(lsqr->vwork_m,lsqr->nwork_m);
     30:   }
     31:   KSPGetVecs(ksp,nw,&lsqr->vwork_m);
     32:   PetscLogObjectParents(ksp,nw,lsqr->vwork_m);

     34:   lsqr->nwork_n = nw = 3;
     35:   if (lsqr->vwork_n) {
     36:     VecDestroyVecs(lsqr->vwork_n,lsqr->nwork_n);
     37:   }
     38:   KSPGetVecs(ksp,nw,&lsqr->vwork_n);
     39:   PetscLogObjectParents(ksp,nw,lsqr->vwork_n);

     41:   return(0);
     42: }

     46: static PetscErrorCode KSPSolve_LSQR(KSP ksp)
     47: {
     49:   PetscInt       i;
     50:   PetscScalar    rho,rhobar,phi,phibar,theta,c,s,tmp,zero = 0.0,mone=-1.0;
     51:   PetscReal      beta,alpha,rnorm;
     52:   Vec            X,B,V,V1,U,U1,TMP,W;
     53:   Mat            Amat,Pmat;
     54:   MatStructure   pflag;
     55:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;
     56:   PetscTruth     diagonalscale;

     59:   PCDiagonalScale(ksp->pc,&diagonalscale);
     60:   if (diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",ksp->type_name);

     62:   PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);

     64:   /* vectors of length m, where system size is mxn */
     65:   B        = ksp->vec_rhs;
     66:   U        = lsqr->vwork_m[0];
     67:   U1       = lsqr->vwork_m[1];

     69:   /* vectors of length n */
     70:   X        = ksp->vec_sol;
     71:   W        = lsqr->vwork_n[0];
     72:   V        = lsqr->vwork_n[1];
     73:   V1       = lsqr->vwork_n[2];

     75:   /* Compute initial residual, temporarily use work vector u */
     76:   if (!ksp->guess_zero) {
     77:     KSP_MatMult(ksp,Amat,X,U);       /*   u <- b - Ax     */
     78:     VecAYPX(U,mone,B);
     79:   } else {
     80:     VecCopy(B,U);            /*   u <- b (x is 0) */
     81:   }

     83:   /* Test for nothing to do */
     84:   VecNorm(U,NORM_2,&rnorm);
     85:   PetscObjectTakeAccess(ksp);
     86:   ksp->its   = 0;
     87:   ksp->rnorm = rnorm;
     88:   PetscObjectGrantAccess(ksp);
     89:   KSPLogResidualHistory(ksp,rnorm);
     90:   KSPMonitor(ksp,0,rnorm);
     91:   (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
     92:   if (ksp->reason) return(0);

     94:   VecCopy(B,U);
     95:   VecNorm(U,NORM_2,&beta);
     96:   tmp = 1.0/beta; VecScale(U,tmp);
     97:   KSP_MatMultTranspose(ksp,Amat,U,V);
     98:   VecNorm(V,NORM_2,&alpha);
     99:   tmp = 1.0/alpha; VecScale(V,tmp);

    101:   VecCopy(V,W);
    102:   VecSet(X,zero);

    104:   phibar = beta;
    105:   rhobar = alpha;
    106:   i = 0;
    107:   do {

    109:     KSP_MatMult(ksp,Amat,V,U1);
    110:     tmp  = -alpha; VecAXPY(U1,tmp,U);
    111:     VecNorm(U1,NORM_2,&beta);
    112:     tmp  = 1.0/beta; VecScale(U1,tmp);

    114:     KSP_MatMultTranspose(ksp,Amat,U1,V1);
    115:     tmp  = -beta; VecAXPY(V1,tmp,V);
    116:     VecNorm(V1,NORM_2,&alpha);
    117:     tmp  = 1.0 / alpha; VecScale(V1,tmp);

    119:     rho    = PetscSqrtScalar(rhobar*rhobar + beta*beta);
    120:     c      = rhobar / rho;
    121:     s      = beta / rho;
    122:     theta  = s * alpha;
    123:     rhobar = - c * alpha;
    124:     phi    = c * phibar;
    125:     phibar = s * phibar;

    127:     tmp  = phi/rho;
    128:     VecAXPY(X,tmp,W);  /*    x <- x + (phi/rho) w   */
    129:     tmp  = -theta/rho;
    130:     VecAYPX(W,tmp,V1); /*    w <- v - (theta/rho) w */

    132:     rnorm = PetscRealPart(phibar);

    134:     PetscObjectTakeAccess(ksp);
    135:     ksp->its++;
    136:     ksp->rnorm = rnorm;
    137:     PetscObjectGrantAccess(ksp);
    138:     KSPLogResidualHistory(ksp,rnorm);
    139:     KSPMonitor(ksp,i+1,rnorm);
    140:     (*ksp->converged)(ksp,i+1,rnorm,&ksp->reason,ksp->cnvP);
    141:     if (ksp->reason) break;
    142:     SWAP(U1,U,TMP);
    143:     SWAP(V1,V,TMP);

    145:     i++;
    146:   } while (i<ksp->max_it);
    147:   if (i >= ksp->max_it && !ksp->reason) {
    148:     ksp->reason = KSP_DIVERGED_ITS;
    149:   }

    151:   /* KSPUnwindPreconditioner(ksp,X,W); */

    153:   return(0);
    154: }

    158: PetscErrorCode KSPDestroy_LSQR(KSP ksp)
    159: {
    160:   KSP_LSQR       *lsqr = (KSP_LSQR*)ksp->data;

    165:   /* Free work vectors */
    166:   if (lsqr->vwork_n) {
    167:     VecDestroyVecs(lsqr->vwork_n,lsqr->nwork_n);
    168:   }
    169:   if (lsqr->vwork_m) {
    170:     VecDestroyVecs(lsqr->vwork_m,lsqr->nwork_m);
    171:   }
    172:   PetscFree(lsqr);
    173:   return(0);
    174: }

    176: /*MC
    177:      KSPLSQR - This implements LSQR

    179:    Options Database Keys:
    180: .   see KSPSolve()

    182:    Level: beginner

    184:    Notes:  This algorithm DOES NOT use a preconditioner. It ignores any preconditioner arguments specified.
    185:            Reference: Paige and Saunders, ACM Transactions on Mathematical Software, Vol 8, pp 43-71, 1982

    187: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP

    189: M*/
    193: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_LSQR(KSP ksp)
    194: {
    195:   KSP_LSQR       *lsqr;

    199:   PetscMalloc(sizeof(KSP_LSQR),&lsqr);
    200:   PetscMemzero(lsqr,sizeof(KSP_LSQR));
    201:   PetscLogObjectMemory(ksp,sizeof(KSP_LSQR));
    202:   ksp->data                      = (void*)lsqr;
    203:   ksp->pc_side                   = PC_LEFT;
    204:   ksp->ops->setup                = KSPSetUp_LSQR;
    205:   ksp->ops->solve                = KSPSolve_LSQR;
    206:   ksp->ops->destroy              = KSPDestroy_LSQR;
    207:   ksp->ops->buildsolution        = KSPDefaultBuildSolution;
    208:   ksp->ops->buildresidual        = KSPDefaultBuildResidual;
    209:   ksp->ops->setfromoptions       = 0;
    210:   ksp->ops->view                 = 0;
    211:   return(0);
    212: }