Actual source code: cgls.c

  1: /*
  2:     This file implements CGLS, the Conjugate Gradient method for Least-Squares problems.
  3: */
  4: #include <petsc/private/kspimpl.h>

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

 12: static PetscErrorCode KSPSetUp_CGLS(KSP ksp)
 13: {
 14:   KSP_CGLS       *cgls = (KSP_CGLS*)ksp->data;

 16:   cgls->nwork_m = 2;
 17:   if (cgls->vwork_m) {
 18:     VecDestroyVecs(cgls->nwork_m,&cgls->vwork_m);
 19:   }

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

 29: static PetscErrorCode KSPSolve_CGLS(KSP ksp)
 30: {
 31:   KSP_CGLS       *cgls = (KSP_CGLS*)ksp->data;
 32:   Mat            A;
 33:   Vec            x,b,r,p,q,ss;
 34:   PetscScalar    beta;
 35:   PetscReal      alpha,gamma,oldgamma;

 37:   KSPGetOperators(ksp,&A,NULL); /* Matrix of the system */

 39:   /* vectors of length n, where system size is mxn */
 40:   x  = ksp->vec_sol; /* Solution vector */
 41:   p  = cgls->vwork_n[0];
 42:   ss  = cgls->vwork_n[1];

 44:   /* vectors of length m, where system size is mxn */
 45:   b  = ksp->vec_rhs; /* Right-hand side vector */
 46:   r  = cgls->vwork_m[0];
 47:   q  = cgls->vwork_m[1];

 49:   /* Minimization with the CGLS method */
 50:   ksp->its = 0;
 51:   ksp->rnorm = 0;
 52:   MatMult(A,x,r);
 53:   VecAYPX(r,-1,b);         /* r_0 = b - A * x_0  */
 54:   KSP_MatMultHermitianTranspose(ksp,A,r,p); /* p_0 = A^T * r_0    */
 55:   VecCopy(p,ss);           /* s_0 = p_0          */
 56:   VecNorm(ss,NORM_2,&gamma);
 57:   KSPCheckNorm(ksp,gamma);
 58:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = gamma;
 59:   (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);
 60:   if (ksp->reason) return 0;
 61:   gamma = gamma*gamma;                          /* gamma = norm2(s)^2 */

 63:   do {
 64:     MatMult(A,p,q);           /* q = A * p               */
 65:     VecNorm(q,NORM_2,&alpha);
 66:     KSPCheckNorm(ksp,alpha);
 67:     alpha = alpha*alpha;                           /* alpha = norm2(q)^2      */
 68:     alpha = gamma / alpha;                         /* alpha = gamma / alpha   */
 69:     VecAXPY(x,alpha,p);       /* x += alpha * p          */
 70:     VecAXPY(r,-alpha,q);      /* r -= alpha * q          */
 71:     KSP_MatMultHermitianTranspose(ksp,A,r,ss); /* ss = A^T * r            */
 72:     oldgamma = gamma;                              /* oldgamma = gamma        */
 73:     VecNorm(ss,NORM_2,&gamma);
 74:     KSPCheckNorm(ksp,gamma);
 75:     ksp->its++;
 76:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = gamma;
 77:     KSPMonitor(ksp,ksp->its,ksp->rnorm);
 78:     (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);
 79:     if (ksp->reason) return 0;
 80:     gamma = gamma*gamma;                           /* gamma = norm2(s)^2      */
 81:     beta = gamma/oldgamma;                         /* beta = gamma / oldgamma */
 82:     VecAYPX(p,beta,ss);       /* p = s + beta * p        */
 83:   } while (ksp->its<ksp->max_it);

 85:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
 86:   return 0;
 87: }

 89: static PetscErrorCode KSPDestroy_CGLS(KSP ksp)
 90: {
 91:   KSP_CGLS       *cgls = (KSP_CGLS*)ksp->data;

 93:   /* Free work vectors */
 94:   if (cgls->vwork_n) {
 95:     VecDestroyVecs(cgls->nwork_n,&cgls->vwork_n);
 96:   }
 97:   if (cgls->vwork_m) {
 98:     VecDestroyVecs(cgls->nwork_m,&cgls->vwork_m);
 99:   }
100:   PetscFree(ksp->data);
101:   return 0;
102: }

104: /*MC
105:      KSPCGLS - Conjugate Gradient method for Least-Squares problems

107:    Level: beginner

109:    Supports non-square (rectangular) matrices.

111:    Notes:
112:     This does not use the preconditioner, so one should probably use KSPLSQR instead.

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

117: M*/
118: PETSC_EXTERN PetscErrorCode KSPCreate_CGLS(KSP ksp)
119: {
120:   KSP_CGLS       *cgls;

122:   PetscNewLog(ksp,&cgls);
123:   ksp->data                = (void*)cgls;
124:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
125:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
126:   ksp->ops->setup          = KSPSetUp_CGLS;
127:   ksp->ops->solve          = KSPSolve_CGLS;
128:   ksp->ops->destroy        = KSPDestroy_CGLS;
129:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
130:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
131:   ksp->ops->setfromoptions = NULL;
132:   ksp->ops->view           = NULL;
133:   return 0;
134: }