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: }