1: /*
  2:     This file implements CGLS, the Conjugate Gradient method for Least-Squares problems.
  3: */
  4: #include <petsc/private/kspimpl.h>        /*I "petscksp.h" I*/
  5: #include <petscksp.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 n */
 11: } KSP_CGLS;

 15: static PetscErrorCode KSPSetUp_CGLS(KSP ksp)
 16: {
 18:   KSP_CGLS       *cgls = (KSP_CGLS*)ksp->data;

 21:   cgls->nwork_m = 2;
 22:   if (cgls->vwork_m) {
 23:     VecDestroyVecs(cgls->nwork_m,&cgls->vwork_m);
 24:   }
 25: 
 26:   cgls->nwork_n = 2;
 27:   if (cgls->vwork_n) {
 28:     VecDestroyVecs(cgls->nwork_n,&cgls->vwork_n);
 29:   }
 30:   KSPCreateVecs(ksp,cgls->nwork_n,&cgls->vwork_n,cgls->nwork_m,&cgls->vwork_m);
 31:   return(0);
 32: }

 36: static PetscErrorCode KSPSolve_CGLS(KSP ksp)
 37: {
 39:   KSP_CGLS       *cgls = (KSP_CGLS*)ksp->data;
 40:   Mat            A;
 41:   Vec            x,b,r,p,q,ss;
 42:   PetscScalar    beta;
 43:   PetscReal      alpha,gamma,oldgamma,rnorm;
 44:   PetscInt       maxiter_ls = 15;
 45: 
 47:   KSPGetOperators(ksp,&A,NULL); /* Matrix of the system */
 48: 
 49:   /* vectors of length n, where system size is mxn */
 50:   x  = ksp->vec_sol; /* Solution vector */
 51:   p  = cgls->vwork_n[0];
 52:   ss  = cgls->vwork_n[1];
 53: 
 54:   /* vectors of length m, where system size is mxn */
 55:   b  = ksp->vec_rhs; /* Right-hand side vector */
 56:   r  = cgls->vwork_m[0];
 57:   q  = cgls->vwork_m[1];
 58: 
 59:   /* Minimization with the CGLS method */
 60:   ksp->its = 0;
 61:   MatMult(A,x,r);
 62:   VecAYPX(r,-1,b);         /* r_0 = b - A * x_0  */
 63:   MatMultTranspose(A,r,p); /* p_0 = A^T * r_0    */
 64:   VecCopy(p,ss);           /* s_0 = p_0          */
 65:   VecNorm(ss,NORM_2,&gamma);
 66:   gamma = gamma*gamma;                          /* gamma = norm2(s)^2 */
 67: 
 68:   do {
 69:     MatMult(A,p,q);           /* q = A * p               */
 70:     VecNorm(q,NORM_2,&alpha);
 71:     alpha = alpha*alpha;                           /* alpha = norm2(q)^2      */
 72:     alpha = gamma / alpha;                         /* alpha = gamma / alpha   */
 73:     VecAXPY(x,alpha,p);       /* x += alpha * p          */
 74:     VecAXPY(r,-alpha,q);      /* r -= alpha * q          */
 75:     MatMultTranspose(A,r,ss); /* ss = A^T * r            */
 76:     oldgamma = gamma;                              /* oldgamma = gamma        */
 77:     VecNorm(ss,NORM_2,&gamma);
 78:     gamma = gamma*gamma;                           /* gamma = norm2(s)^2      */
 79:     beta = gamma/oldgamma;                         /* beta = gamma / oldgamma */
 80:     VecAYPX(p,beta,ss);       /* p = s + beta * p        */
 81:     ksp->its ++;
 82:     ksp->rnorm = gamma;
 83:   } while (ksp->its<ksp->max_it);
 84: 
 85:   MatMult(A,x,r);
 86:   VecAXPY(r,-1,b);
 87:   VecNorm(r,NORM_2,&rnorm);
 88:   ksp->rnorm = rnorm;
 89:   KSPMonitor(ksp,ksp->its,rnorm);
 90:   (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
 91:   if (ksp->its>=maxiter_ls && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
 92:   return(0);
 93: }

 97: static PetscErrorCode KSPDestroy_CGLS(KSP ksp)
 98: {
 99:   KSP_CGLS       *cgls = (KSP_CGLS*)ksp->data;

103:   /* Free work vectors */
104:   if (cgls->vwork_n) {
105:     VecDestroyVecs(cgls->nwork_n,&cgls->vwork_n);
106:   }
107:   if (cgls->vwork_m) {
108:     VecDestroyVecs(cgls->nwork_m,&cgls->vwork_m);
109:   }
110:   PetscFree(ksp->data);
111:   return(0);
112: }

114: /*MC
115:      KSCGLS - Conjugate Gradient method for Least-Squares problems

117:    Level: beginner

119:    Supports non-square (rectangular) matrices.

121: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
122:            KSPCGSetType(), KSPCGUseSingleReduction(), KSPPIPECG, KSPGROPPCG

124: M*/

128: PETSC_EXTERN PetscErrorCode KSPCreate_CGLS(KSP ksp)
129: {
131:   KSP_CGLS       *cgls;
132: 
134:   PetscNewLog(ksp,&cgls);
135:   ksp->data                = (void*)cgls;
136:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
137:   ksp->ops->setup          = KSPSetUp_CGLS;
138:   ksp->ops->solve          = KSPSolve_CGLS;
139:   ksp->ops->destroy        = KSPDestroy_CGLS;
140:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
141:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
142:   ksp->ops->setfromoptions = 0;
143:   ksp->ops->view           = 0;
144: #if defined(PETSC_USE_COMPLEX)
145:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"This is not supported for complex numbers");
146: #endif
147:   return(0);
148: }