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) VecDestroyVecs(cgls->nwork_m, &cgls->vwork_m);

 19:   cgls->nwork_n = 2;
 20:   if (cgls->vwork_n) VecDestroyVecs(cgls->nwork_n, &cgls->vwork_n);
 21:   KSPCreateVecs(ksp, cgls->nwork_n, &cgls->vwork_n, cgls->nwork_m, &cgls->vwork_m);
 22:   return 0;
 23: }

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

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

 35:   /* vectors of length n, where system size is mxn */
 36:   x  = ksp->vec_sol; /* Solution vector */
 37:   p  = cgls->vwork_n[0];
 38:   ss = cgls->vwork_n[1];

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

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

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

 81:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
 82:   return 0;
 83: }

 85: static PetscErrorCode KSPDestroy_CGLS(KSP ksp)
 86: {
 87:   KSP_CGLS *cgls = (KSP_CGLS *)ksp->data;

 89:   /* Free work vectors */
 90:   if (cgls->vwork_n) VecDestroyVecs(cgls->nwork_n, &cgls->vwork_n);
 91:   if (cgls->vwork_m) VecDestroyVecs(cgls->nwork_m, &cgls->vwork_m);
 92:   PetscFree(ksp->data);
 93:   return 0;
 94: }

 96: /*MC
 97:      KSPCGLS - Conjugate Gradient method for Least-Squares problems. Supports non-square (rectangular) matrices.

 99:    Level: beginner

101:    Note:
102:    This does not use the preconditioner, so one should probably use `KSPLSQR` instead.

104: .seealso: [](chapter_ksp), `KSPLSQR`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
105:           `KSPCGSetType()`, `KSPCGUseSingleReduction()`, `KSPPIPECG`, `KSPGROPPCG`
106: M*/
107: PETSC_EXTERN PetscErrorCode KSPCreate_CGLS(KSP ksp)
108: {
109:   KSP_CGLS *cgls;

111:   PetscNew(&cgls);
112:   ksp->data = (void *)cgls;
113:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3);
114:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
115:   ksp->ops->setup          = KSPSetUp_CGLS;
116:   ksp->ops->solve          = KSPSolve_CGLS;
117:   ksp->ops->destroy        = KSPDestroy_CGLS;
118:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
119:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
120:   ksp->ops->setfromoptions = NULL;
121:   ksp->ops->view           = NULL;
122:   return 0;
123: }