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:   PetscFunctionBegin;
 17:   cgls->nwork_m = 2;
 18:   if (cgls->vwork_m) PetscCall(VecDestroyVecs(cgls->nwork_m, &cgls->vwork_m));

 20:   cgls->nwork_n = 2;
 21:   if (cgls->vwork_n) PetscCall(VecDestroyVecs(cgls->nwork_n, &cgls->vwork_n));
 22:   PetscCall(KSPCreateVecs(ksp, cgls->nwork_n, &cgls->vwork_n, cgls->nwork_m, &cgls->vwork_m));
 23:   PetscFunctionReturn(PETSC_SUCCESS);
 24: }

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

 34:   PetscFunctionBegin;
 35:   PetscCall(KSPGetOperators(ksp, &A, NULL)); /* Matrix of the system */

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

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

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

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

 83:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
 84:   PetscFunctionReturn(PETSC_SUCCESS);
 85: }

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

 91:   PetscFunctionBegin;
 92:   /* Free work vectors */
 93:   if (cgls->vwork_n) PetscCall(VecDestroyVecs(cgls->nwork_n, &cgls->vwork_n));
 94:   if (cgls->vwork_m) PetscCall(VecDestroyVecs(cgls->nwork_m, &cgls->vwork_m));
 95:   PetscCall(PetscFree(ksp->data));
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

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

102:    Level: beginner

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

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

114:   PetscFunctionBegin;
115:   PetscCall(PetscNew(&cgls));
116:   ksp->data = (void *)cgls;
117:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
118:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
119:   ksp->ops->setup          = KSPSetUp_CGLS;
120:   ksp->ops->solve          = KSPSolve_CGLS;
121:   ksp->ops->destroy        = KSPDestroy_CGLS;
122:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
123:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
124:   ksp->ops->setfromoptions = NULL;
125:   ksp->ops->view           = NULL;
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }