Actual source code: tsirm.c

      1: #include <petsc/private/kspimpl.h>

      3: typedef struct {
      4:   PetscReal tol_ls;
      5:   PetscInt  size_ls, maxiter_ls, cgls, size, Istart, Iend;
      6:   Mat       A, S;
      7:   Vec       Alpha, r;
      8: } KSP_TSIRM;

     10: static PetscErrorCode KSPSetUp_TSIRM(KSP ksp)
     11: {
     12:   KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;

     14:   PetscFunctionBegin;
     15:   /* Matrix of the system */
     16:   PetscCall(KSPGetOperators(ksp, &tsirm->A, NULL));    /* Matrix of the system   */
     17:   PetscCall(MatGetSize(tsirm->A, &tsirm->size, NULL)); /* Size of the system     */
     18:   PetscCall(MatGetOwnershipRange(tsirm->A, &tsirm->Istart, &tsirm->Iend));

     20:   /* Matrix S of residuals */
     21:   PetscCall(MatCreate(PetscObjectComm((PetscObject)ksp), &tsirm->S));
     22:   PetscCall(MatSetSizes(tsirm->S, tsirm->Iend - tsirm->Istart, PETSC_DECIDE, tsirm->size, tsirm->size_ls));
     23:   PetscCall(MatSetType(tsirm->S, MATDENSE));
     24:   PetscCall(MatSetUp(tsirm->S));

     26:   /* Residual and vector Alpha computed in the minimization step */
     27:   PetscCall(MatCreateVecs(tsirm->S, &tsirm->Alpha, &tsirm->r));
     28:   PetscFunctionReturn(PETSC_SUCCESS);
     29: }

     31: static PetscErrorCode KSPSolve_TSIRM(KSP ksp)
     32: {
     33:   KSP_TSIRM   *tsirm = (KSP_TSIRM *)ksp->data;
     34:   KSP          sub_ksp;
     35:   PC           pc;
     36:   Mat          AS = NULL;
     37:   Vec          x, b;
     38:   PetscScalar *array;
     39:   PetscReal    norm = 20;
     40:   PetscInt     i, *ind_row, first_iteration = 1, its = 0, total = 0, col = 0;
     41:   KSP          ksp_min; /* KSP for minimization */
     42:   PC           pc_min;  /* PC for minimization */
     43:   PetscBool    isksp;

     45:   PetscFunctionBegin;
     46:   x = ksp->vec_sol; /* Solution vector        */
     47:   b = ksp->vec_rhs; /* Right-hand side vector */

     49:   /* Row indexes (these indexes are global) */
     50:   PetscCall(PetscMalloc1(tsirm->Iend - tsirm->Istart, &ind_row));
     51:   for (i = 0; i < tsirm->Iend - tsirm->Istart; i++) ind_row[i] = i + tsirm->Istart;

     53:   /* Inner solver */
     54:   PetscCall(KSPGetPC(ksp, &pc));
     55:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCKSP, &isksp));
     56:   PetscCheck(isksp, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "PC must be of type PCKSP");
     57:   PetscCall(PCKSPGetKSP(pc, &sub_ksp));

     59:   /* previously it seemed good but with SNES it seems not good... */
     60:   PetscCall(KSP_MatMult(sub_ksp, tsirm->A, x, tsirm->r));
     61:   PetscCall(VecAXPY(tsirm->r, -1, b));
     62:   PetscCall(VecNorm(tsirm->r, NORM_2, &norm));
     63:   KSPCheckNorm(ksp, norm);
     64:   ksp->its = 0;
     65:   PetscCall(KSPConvergedDefault(ksp, ksp->its, norm, &ksp->reason, ksp->cnvP));
     66:   PetscCall(KSPMonitor(ksp, ksp->its, norm));
     67:   PetscCall(KSPSetInitialGuessNonzero(sub_ksp, PETSC_TRUE));
     68:   do {
     69:     for (col = 0; col < tsirm->size_ls && ksp->reason == KSP_CONVERGED_ITERATING; col++) {
     70:       /* Solve (inner iteration) */
     71:       PetscCall(KSPSolve(sub_ksp, b, x));
     72:       PetscCall(KSPGetIterationNumber(sub_ksp, &its));
     73:       total += its;

     75:       /* Build S^T */
     76:       PetscCall(VecGetArray(x, &array));
     77:       PetscCall(MatSetValues(tsirm->S, tsirm->Iend - tsirm->Istart, ind_row, 1, &col, array, INSERT_VALUES));
     78:       PetscCall(VecRestoreArray(x, &array));

     80:       PetscCall(KSPGetResidualNorm(sub_ksp, &norm));
     81:       ksp->rnorm = norm;
     82:       ksp->its++;
     83:       PetscCall(KSPConvergedDefault(ksp, ksp->its, norm, &ksp->reason, ksp->cnvP));
     84:       PetscCall(KSPMonitor(ksp, ksp->its, norm));
     85:     }

     87:     /* Minimization step */
     88:     if (ksp->reason == KSP_CONVERGED_ITERATING) {
     89:       PetscCall(MatAssemblyBegin(tsirm->S, MAT_FINAL_ASSEMBLY));
     90:       PetscCall(MatAssemblyEnd(tsirm->S, MAT_FINAL_ASSEMBLY));
     91:       if (first_iteration) {
     92:         PetscCall(MatMatMult(tsirm->A, tsirm->S, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AS));
     93:         first_iteration = 0;
     94:       } else {
     95:         PetscCall(MatMatMult(tsirm->A, tsirm->S, MAT_REUSE_MATRIX, PETSC_DETERMINE, &AS));
     96:       }

     98:       /* CGLS or LSQR method to minimize the residuals*/
     99:       PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &ksp_min));
    100:       if (tsirm->cgls) {
    101:         PetscCall(KSPSetType(ksp_min, KSPCGLS));
    102:       } else {
    103:         PetscCall(KSPSetType(ksp_min, KSPLSQR));
    104:       }
    105:       PetscCall(KSPSetOperators(ksp_min, AS, AS));
    106:       PetscCall(KSPSetTolerances(ksp_min, tsirm->tol_ls, PETSC_CURRENT, PETSC_CURRENT, tsirm->maxiter_ls));
    107:       PetscCall(KSPGetPC(ksp_min, &pc_min));
    108:       PetscCall(PCSetType(pc_min, PCNONE));
    109:       PetscCall(KSPSolve(ksp_min, b, tsirm->Alpha)); /* Find Alpha such that ||AS Alpha = b|| */
    110:       PetscCall(KSPDestroy(&ksp_min));
    111:       /* Apply minimization */
    112:       PetscCall(MatMult(tsirm->S, tsirm->Alpha, x)); /* x = S * Alpha */
    113:     }
    114:   } while (ksp->its < ksp->max_it && !ksp->reason);
    115:   PetscCall(MatDestroy(&AS));
    116:   PetscCall(PetscFree(ind_row));
    117:   ksp->its = total;
    118:   PetscFunctionReturn(PETSC_SUCCESS);
    119: }

    121: static PetscErrorCode KSPSetFromOptions_TSIRM(KSP ksp, PetscOptionItems PetscOptionsObject)
    122: {
    123:   KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;

    125:   PetscFunctionBegin;
    126:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP TSIRM options");
    127:   PetscCall(PetscOptionsInt("-ksp_tsirm_cgls", "Method used for the minimization step", "", tsirm->cgls, &tsirm->cgls, NULL)); /*0:LSQR, 1:CGLS*/
    128:   PetscCall(PetscOptionsReal("-ksp_tsirm_tol_ls", "Tolerance threshold for the minimization step", "", tsirm->tol_ls, &tsirm->tol_ls, NULL));
    129:   PetscCall(PetscOptionsInt("-ksp_tsirm_max_it_ls", "Maximum number of iterations for the minimization step", "", tsirm->maxiter_ls, &tsirm->maxiter_ls, NULL));
    130:   PetscCall(PetscOptionsInt("-ksp_tsirm_size_ls", "Number of residuals for minimization", "", tsirm->size_ls, &tsirm->size_ls, NULL));
    131:   PetscOptionsHeadEnd();
    132:   PetscFunctionReturn(PETSC_SUCCESS);
    133: }

    135: static PetscErrorCode KSPDestroy_TSIRM(KSP ksp)
    136: {
    137:   KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;

    139:   PetscFunctionBegin;
    140:   PetscCall(MatDestroy(&tsirm->S));
    141:   PetscCall(VecDestroy(&tsirm->Alpha));
    142:   PetscCall(VecDestroy(&tsirm->r));
    143:   PetscCall(PetscFree(ksp->data));
    144:   PetscFunctionReturn(PETSC_SUCCESS);
    145: }

    147: /*MC
    148:    KSPTSIRM - Implements the two-stage iteration with least-squares residual minimization method {cite}`couturier2016tsirm`

    150:    Options Database Keys:
    151: +  -ksp_ksp_type <solver>        - the type of the inner solver (GMRES or any of its variants for instance)
    152: .  -ksp_pc_type <preconditioner> - the type of the preconditioner applied to the inner solver
    153: .  -ksp_ksp_max_it <maxits>      - the maximum number of inner iterations (iterations of the inner solver)
    154: .  -ksp_ksp_rtol <tol>           - sets the relative convergence tolerance of the inner solver
    155: .  -ksp_tsirm_cgls <number>      - if 1 use CGLS solver in the minimization step, otherwise use LSQR solver
    156: .  -ksp_tsirm_max_it_ls <maxits> - the maximum number of iterations for the least-squares minimization solver
    157: .  -ksp_tsirm_tol_ls <tol>       - sets the convergence tolerance of the least-squares minimization solver
    158: -  -ksp_tsirm_size_ls <size>     - the number of residuals for the least-squares minimization step

    160:    Level: advanced

    162:    Notes:
    163:    `KSPTSIRM` is a two-stage iteration method for solving large sparse linear systems of the form $Ax=b$. The main idea behind this new
    164:    method is the use a least-squares residual minimization to improve the convergence of Krylov based iterative methods, typically those of GMRES variants.
    165:    The principle of `TSIRM` algorithm  is to build an outer iteration over a Krylov method, called the inner solver, and to frequently store the current residual
    166:    computed by the given Krylov method in a matrix of residuals S. After a few outer iterations, a least-squares minimization step is applied on the matrix
    167:    composed by the saved residuals, in order to compute a better solution and to make new iterations if required.
    168:    The minimization step consists in solving the least-squares problem $\min||b-ASa||$ to find 'a' which minimizes the
    169:    residuals $(b-AS)$. The minimization step is performed using two solvers of linear least-squares problems: `KSPCGLS` or `KSPLSQR`. A new solution x with
    170:    a minimal residual is computed with $x=Sa$.

    172:    Defaults to 30 iterations for the inner solve, use option `-ksp_ksp_max_it <it>` to change it.

    174:    Contributed by:
    175:    Lilia Ziane Khodja

    177: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPLGMRES`,
    178:           `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
    179:           `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
    180:           `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPSetPCSide()`
    181: M*/
    182: PETSC_EXTERN PetscErrorCode KSPCreate_TSIRM(KSP ksp)
    183: {
    184:   KSP_TSIRM *tsirm;
    185:   PC         pc;
    186:   KSP        sub_ksp;

    188:   PetscFunctionBegin;
    189:   PetscCall(PetscNew(&tsirm));
    190:   ksp->data = (void *)tsirm;
    191: #if defined(PETSC_USE_REAL_SINGLE)
    192:   tsirm->tol_ls = 1e-25;
    193: #else
    194:   tsirm->tol_ls = 1e-50;
    195: #endif
    196:   tsirm->size_ls    = 12;
    197:   tsirm->maxiter_ls = 15;
    198:   tsirm->cgls       = 0;
    199:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
    200:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1));
    201:   ksp->ops->setup          = KSPSetUp_TSIRM;
    202:   ksp->ops->solve          = KSPSolve_TSIRM;
    203:   ksp->ops->destroy        = KSPDestroy_TSIRM;
    204:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    205:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    206:   ksp->ops->setfromoptions = KSPSetFromOptions_TSIRM;
    207:   ksp->ops->view           = NULL;

    209:   PetscCall(KSPGetPC(ksp, &pc));
    210:   PetscCall(PCSetType(pc, PCKSP));
    211:   PetscCall(PCKSPGetKSP(pc, &sub_ksp));
    212:   PetscCall(KSPSetTolerances(sub_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 30));
    213: #if defined(PETSC_USE_COMPLEX)
    214:   SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "This is not supported for complex numbers");
    215: #else
    216:   PetscFunctionReturn(PETSC_SUCCESS);
    217: #endif
    218: }