Actual source code: cgne.c

      1: /*
      2:     cgimpl.h defines the simple data structured used to store information
      3:     related to the type of matrix (e.g. complex symmetric) being solved and
      4:     data used during the optional Lanczos process used to compute eigenvalues
      5: */
      6: #include <../src/ksp/ksp/impls/cg/cgimpl.h>
      7: extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP, PetscReal *, PetscReal *);
      8: extern PetscErrorCode KSPComputeEigenvalues_CG(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);

     10: static PetscErrorCode KSPCGSetType_CGNE(KSP ksp, KSPCGType type)
     11: {
     12:   KSP_CG *cg = (KSP_CG *)ksp->data;

     14:   PetscFunctionBegin;
     15:   cg->type = type;
     16:   PetscFunctionReturn(PETSC_SUCCESS);
     17: }

     19: static PetscErrorCode KSPSetUp_CGNE(KSP ksp)
     20: {
     21:   KSP_CG  *cgP   = (KSP_CG *)ksp->data;
     22:   PetscInt maxit = ksp->max_it;

     24:   PetscFunctionBegin;
     25:   /* get work vectors needed by CGNE */
     26:   PetscCall(KSPSetWorkVecs(ksp, 4));

     28:   /*
     29:      If user requested computations of eigenvalues then allocate work space needed
     30:   */
     31:   if (ksp->calc_sings) {
     32:     /* get space to store tridiagonal matrix for Lanczos */
     33:     PetscCall(PetscMalloc4(maxit, &cgP->e, maxit, &cgP->d, maxit, &cgP->ee, maxit, &cgP->dd));

     35:     ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
     36:     ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_CG;
     37:   }
     38:   PetscFunctionReturn(PETSC_SUCCESS);
     39: }

     41: static PetscErrorCode KSPSolve_CGNE(KSP ksp)
     42: {
     43:   PetscInt    i, stored_max_it, eigs;
     44:   PetscScalar dpi, a = 1.0, beta, betaold = 1.0, b = 0, *e = NULL, *d = NULL;
     45:   PetscReal   dp = 0.0;
     46:   Vec         X, B, Z, R, P, T;
     47:   KSP_CG     *cg;
     48:   Mat         Amat, Pmat;
     49:   PetscBool   diagonalscale, transpose_pc;

     51:   PetscFunctionBegin;
     52:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
     53:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
     54:   PetscCall(PCApplyTransposeExists(ksp->pc, &transpose_pc));

     56:   cg            = (KSP_CG *)ksp->data;
     57:   eigs          = ksp->calc_sings;
     58:   stored_max_it = ksp->max_it;
     59:   X             = ksp->vec_sol;
     60:   B             = ksp->vec_rhs;
     61:   R             = ksp->work[0];
     62:   Z             = ksp->work[1];
     63:   P             = ksp->work[2];
     64:   T             = ksp->work[3];

     66: #define VecXDot(x, y, a) (cg->type == KSP_CG_HERMITIAN ? VecDot(x, y, a) : VecTDot(x, y, a))

     68:   if (eigs) {
     69:     e    = cg->e;
     70:     d    = cg->d;
     71:     e[0] = 0.0;
     72:   }
     73:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));

     75:   ksp->its = 0;
     76:   PetscCall(KSP_MatMultTranspose(ksp, Amat, B, T));
     77:   if (!ksp->guess_zero) {
     78:     PetscCall(KSP_MatMult(ksp, Amat, X, P));
     79:     PetscCall(KSP_MatMultTranspose(ksp, Amat, P, R));
     80:     PetscCall(VecAYPX(R, -1.0, T));
     81:   } else {
     82:     PetscCall(VecCopy(T, R)); /*     r <- b (x is 0) */
     83:   }
     84:   if (transpose_pc) {
     85:     PetscCall(KSP_PCApplyTranspose(ksp, R, T));
     86:   } else {
     87:     PetscCall(KSP_PCApply(ksp, R, T));
     88:   }
     89:   PetscCall(KSP_PCApply(ksp, T, Z));

     91:   if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
     92:     PetscCall(VecNorm(Z, NORM_2, &dp)); /*    dp <- z'*z       */
     93:   } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
     94:     PetscCall(VecNorm(R, NORM_2, &dp)); /*    dp <- r'*r       */
     95:   } else if (ksp->normtype == KSP_NORM_NATURAL) {
     96:     PetscCall(VecXDot(Z, R, &beta));
     97:     KSPCheckDot(ksp, beta);
     98:     dp = PetscSqrtReal(PetscAbsScalar(beta));
     99:   } else dp = 0.0;
    100:   PetscCall(KSPLogResidualHistory(ksp, dp));
    101:   PetscCall(KSPMonitor(ksp, 0, dp));
    102:   ksp->rnorm = dp;
    103:   PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); /* test for convergence */
    104:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

    106:   i = 0;
    107:   do {
    108:     ksp->its = i + 1;
    109:     PetscCall(VecXDot(Z, R, &beta)); /*     beta <- r'z     */
    110:     KSPCheckDot(ksp, beta);
    111:     if (beta == 0.0) {
    112:       ksp->reason = KSP_CONVERGED_ATOL;
    113:       PetscCall(PetscInfo(ksp, "converged due to beta = 0\n"));
    114:       break;
    115: #if !defined(PETSC_USE_COMPLEX)
    116:     } else if (beta < 0.0) {
    117:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
    118:       PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
    119:       break;
    120: #endif
    121:     }
    122:     if (!i) {
    123:       PetscCall(VecCopy(Z, P)); /*     p <- z          */
    124:       b = 0.0;
    125:     } else {
    126:       b = beta / betaold;
    127:       if (eigs) {
    128:         PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
    129:         e[i] = PetscSqrtReal(PetscAbsScalar(b)) / a;
    130:       }
    131:       PetscCall(VecAYPX(P, b, Z)); /*     p <- z + b* p   */
    132:     }
    133:     betaold = beta;
    134:     PetscCall(KSP_MatMult(ksp, Amat, P, T));
    135:     PetscCall(KSP_MatMultTranspose(ksp, Amat, T, Z));
    136:     PetscCall(VecXDot(P, Z, &dpi)); /*     dpi <- z'p      */
    137:     KSPCheckDot(ksp, dpi);
    138:     a = beta / dpi; /*     a = beta/p'z    */
    139:     if (eigs) d[i] = PetscSqrtReal(PetscAbsScalar(b)) * e[i] + 1.0 / a;
    140:     PetscCall(VecAXPY(X, a, P));  /*     x <- x + ap     */
    141:     PetscCall(VecAXPY(R, -a, Z)); /*     r <- r - az     */
    142:     if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
    143:       if (transpose_pc) {
    144:         PetscCall(KSP_PCApplyTranspose(ksp, R, T));
    145:       } else {
    146:         PetscCall(KSP_PCApply(ksp, R, T));
    147:       }
    148:       PetscCall(KSP_PCApply(ksp, T, Z));
    149:       PetscCall(VecNorm(Z, NORM_2, &dp)); /*    dp <- z'*z       */
    150:     } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
    151:       PetscCall(VecNorm(R, NORM_2, &dp));
    152:     } else if (ksp->normtype == KSP_NORM_NATURAL) {
    153:       dp = PetscSqrtReal(PetscAbsScalar(beta));
    154:     } else dp = 0.0;
    155:     ksp->rnorm = dp;
    156:     PetscCall(KSPLogResidualHistory(ksp, dp));
    157:     PetscCall(KSPMonitor(ksp, i + 1, dp));
    158:     PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
    159:     if (ksp->reason) break;
    160:     if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
    161:       if (transpose_pc) {
    162:         PetscCall(KSP_PCApplyTranspose(ksp, R, T));
    163:       } else {
    164:         PetscCall(KSP_PCApply(ksp, R, T));
    165:       }
    166:       PetscCall(KSP_PCApply(ksp, T, Z));
    167:     }
    168:     i++;
    169:   } while (i < ksp->max_it);
    170:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    171:   PetscFunctionReturn(PETSC_SUCCESS);
    172: }

    174: /*MC
    175:    KSPCGNE - Applies the preconditioned conjugate gradient method to the normal equations
    176:              without explicitly forming $A^T A$.

    178:    Options Database Key:
    179: .   -ksp_cg_type <Hermitian or symmetric - (for complex matrices only) indicates the matrix is Hermitian or symmetric

    181:    Level: beginner

    183:    Notes:
    184:    Eigenvalue computation routines including `KSPSetComputeEigenvalues()` and `KSPComputeEigenvalues()` will return information about the
    185:    spectrum of $A^T*A$, rather than $A$.

    187:    `KSPCGNE` is a general-purpose non-symmetric method. It works well when the singular values are much better behaved than
    188:    eigenvalues. A unitary matrix is a classic example where `KSPCGNE` converges in one iteration, but `KSPGMRES` and `KSPCGS` need N
    189:    iterations, see {cite}`nachtigal90`. If you intend to solve least squares problems, use `KSPLSQR`.

    191:    This is NOT a different algorithm than used with `KSPCG`, it merely uses that algorithm with the
    192:    matrix defined by $A^T A$ and preconditioner defined by $B^T B$ where $B$ is the preconditioner for $A$.

    194:    See `PETSCREGRESSORLINEAR` for the PETSc toolkit for solving linear regression problems including least squares.

    196:    This method requires that one be able to apply the transpose of the preconditioner and operator
    197:    as well as the operator and preconditioner. If the transpose of the preconditioner is not available then
    198:    the preconditioner is used in its place so one ends up preconditioning $A^T A$ with $B B$. Seems odd?

    200:    This only supports left preconditioning.

    202:    Developer Note:
    203:    This object is subclassed off of `KSPCG`, see the source code in src/ksp/ksp/impls/cg for comments on the structure of the code

    205: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPLSQR`, `KSPCGLS`,
    206:           `KSPCGSetType()`, `KSPBICG`, `KSPSetComputeEigenvalues()`, `KSPComputeEigenvalues()`, `PETSCREGRESSORLINEAR`
    207: M*/

    209: PETSC_EXTERN PetscErrorCode KSPCreate_CGNE(KSP ksp)
    210: {
    211:   KSP_CG *cg;

    213:   PetscFunctionBegin;
    214:   PetscCall(PetscNew(&cg));
    215: #if !defined(PETSC_USE_COMPLEX)
    216:   cg->type = KSP_CG_SYMMETRIC;
    217: #else
    218:   cg->type = KSP_CG_HERMITIAN;
    219: #endif
    220:   ksp->data = (void *)cg;
    221:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
    222:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
    223:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
    224:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

    226:   /*
    227:        Sets the functions that are associated with this data structure
    228:        (in C++ this is the same as defining virtual functions)
    229:   */
    230:   ksp->ops->setup          = KSPSetUp_CGNE;
    231:   ksp->ops->solve          = KSPSolve_CGNE;
    232:   ksp->ops->destroy        = KSPDestroy_CG;
    233:   ksp->ops->view           = KSPView_CG;
    234:   ksp->ops->setfromoptions = KSPSetFromOptions_CG;
    235:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    236:   ksp->ops->buildresidual  = KSPBuildResidualDefault;

    238:   /*
    239:       Attach the function KSPCGSetType_CGNE() to this object. The routine
    240:       KSPCGSetType() checks for this attached function and calls it if it finds
    241:       it. (Sort of like a dynamic member function that can be added at run time
    242:   */
    243:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetType_C", KSPCGSetType_CGNE));
    244:   PetscFunctionReturn(PETSC_SUCCESS);
    245: }