Actual source code: tsirm.c

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

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

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

     18:   /* Initialization */
     19: #if defined(PETSC_USE_REAL_SINGLE)
     20:   tsirm->tol_ls     = 1e-25;
     21: #else
     22:   tsirm->tol_ls     = 1e-50;
     23: #endif
     24:   tsirm->size_ls    = 12;
     25:   tsirm->maxiter_ls = 15;
     26:   tsirm->cgls       = 0;

     28:   /* Matrix of the system */
     29:   KSPGetOperators(ksp,&tsirm->A,NULL);    /* Matrix of the system   */
     30:   MatGetSize(tsirm->A,&tsirm->size,NULL); /* Size of the system     */
     31:   MatGetOwnershipRange(tsirm->A,&tsirm->Istart,&tsirm->Iend);

     33:   /* Matrix S of residuals */
     34:   MatCreate(PETSC_COMM_WORLD,&tsirm->S);
     35:   MatSetSizes(tsirm->S,tsirm->Iend-tsirm->Istart,PETSC_DECIDE,tsirm->size,tsirm->size_ls);
     36:   MatSetType(tsirm->S,MATDENSE);
     37:   MatSetUp(tsirm->S);

     39:   /* Residual and vector Alpha computed in the minimization step */
     40:   MatCreateVecs(tsirm->S,&tsirm->Alpha,&tsirm->r);
     41:   return(0);
     42: }

     44: PetscErrorCode KSPSolve_TSIRM(KSP ksp)
     45: {
     47:   KSP_TSIRM      *tsirm = (KSP_TSIRM*)ksp->data;
     48:   KSP            sub_ksp;
     49:   PC             pc;
     50:   Mat            AS = NULL;
     51:   Vec            x,b;
     52:   PetscScalar    *array;
     53:   PetscReal      norm = 20;
     54:   PetscInt       i,*ind_row,first_iteration = 1,its = 0,total = 0,col = 0;
     55:   PetscInt       restart = 30;
     56:   KSP            ksp_min;  /* KSP for minimization */
     57:   PC             pc_min;    /* PC for minimization */
     58:   PetscBool      isksp;

     61:   x = ksp->vec_sol; /* Solution vector        */
     62:   b = ksp->vec_rhs; /* Right-hand side vector */

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

     68:   /* Inner solver */
     69:   KSPGetPC(ksp,&pc);
     70:   PetscObjectTypeCompare((PetscObject)pc,PCKSP,&isksp);
     71:   if (!isksp) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"PC must be of type PCKSP");
     72:   PCKSPGetKSP(pc,&sub_ksp);
     73:   KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,restart);

     75:   /* previously it seemed good but with SNES it seems not good... */
     76:   KSP_MatMult(sub_ksp,tsirm->A,x,tsirm->r);
     77:   VecAXPY(tsirm->r,-1,b);
     78:   VecNorm(tsirm->r,NORM_2,&norm);
     79:   KSPCheckNorm(ksp,norm);
     80:   ksp->its = 0;
     81:   KSPConvergedDefault(ksp,ksp->its,norm,&ksp->reason,ksp->cnvP);
     82:   KSPSetInitialGuessNonzero(sub_ksp,PETSC_TRUE);
     83:   do {
     84:     for (col=0;col<tsirm->size_ls && ksp->reason==0;col++) {
     85:       /* Solve (inner iteration) */
     86:       KSPSolve(sub_ksp,b,x);
     87:       KSPGetIterationNumber(sub_ksp,&its);
     88:       total += its;

     90:       /* Build S^T */
     91:       VecGetArray(x,&array);
     92:       MatSetValues(tsirm->S,tsirm->Iend-tsirm->Istart,ind_row,1,&col,array,INSERT_VALUES);
     93:       VecRestoreArray(x,&array);

     95:       KSPGetResidualNorm(sub_ksp,&norm);
     96:       ksp->rnorm = norm;
     97:       ksp->its ++;
     98:       KSPConvergedDefault(ksp,ksp->its,norm,&ksp->reason,ksp->cnvP);
     99:       KSPMonitor(ksp,ksp->its,norm);
    100:     }

    102:     /* Minimization step */
    103:     if (!ksp->reason) {
    104:       MatAssemblyBegin(tsirm->S,MAT_FINAL_ASSEMBLY);
    105:       MatAssemblyEnd(tsirm->S,MAT_FINAL_ASSEMBLY);
    106:       if (first_iteration) {
    107:         MatMatMult(tsirm->A,tsirm->S,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
    108:         first_iteration = 0;
    109:       } else {
    110:         MatMatMult(tsirm->A,tsirm->S,MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
    111:       }

    113:       /* CGLS or LSQR method to minimize the residuals*/

    115:       KSPCreate(PETSC_COMM_WORLD,&ksp_min);
    116:       if (tsirm->cgls) {
    117:         KSPSetType(ksp_min,KSPCGLS);
    118:       } else {
    119:         KSPSetType(ksp_min,KSPLSQR);
    120:       }
    121:       KSPSetOperators(ksp_min,AS,AS);
    122:       KSPSetTolerances(ksp_min,tsirm->tol_ls,PETSC_DEFAULT,PETSC_DEFAULT,tsirm->maxiter_ls);
    123:       KSPGetPC(ksp_min,&pc_min);
    124:       PCSetType(pc_min,PCNONE);
    125:       KSPSolve(ksp_min,b,tsirm->Alpha);    /* Find Alpha such that ||AS Alpha = b|| */
    126:       KSPDestroy(&ksp_min);
    127:       /* Apply minimization */
    128:       MatMult(tsirm->S,tsirm->Alpha,x); /* x = S * Alpha */
    129:     }
    130:   } while (ksp->its<ksp->max_it && !ksp->reason);
    131:   MatDestroy(&AS);
    132:   PetscFree(ind_row);
    133:   ksp->its = total;
    134:   return(0);
    135: }

    137: PetscErrorCode KSPSetFromOptions_TSIRM(PetscOptionItems *PetscOptionsObject,KSP ksp)
    138: {
    140:   KSP_TSIRM      *tsirm = (KSP_TSIRM*)ksp->data;

    143:   PetscOptionsHead(PetscOptionsObject,"KSP TSIRM options");
    144:   PetscOptionsInt("-ksp_tsirm_cgls","Method used for the minimization step","",tsirm->cgls,&tsirm->cgls,NULL); /*0:LSQR, 1:CGLS*/
    145:   PetscOptionsReal("-ksp_tsirm_tol_ls","Tolerance threshold for the minimization step","",tsirm->tol_ls,&tsirm->tol_ls,NULL);
    146:   PetscOptionsInt("-ksp_tsirm_max_it_ls","Maximum number of iterations for the minimization step","",tsirm->maxiter_ls,&tsirm->maxiter_ls,NULL);
    147:   PetscOptionsInt("-ksp_tsirm_size_ls","Number of residuals for minimization","",tsirm->size_ls,&tsirm->size_ls,NULL);
    148:   PetscOptionsTail();
    149:   return(0);
    150: }

    152: PetscErrorCode KSPDestroy_TSIRM(KSP ksp)
    153: {
    154:   KSP_TSIRM       *tsirm = (KSP_TSIRM*)ksp->data;

    158:   MatDestroy(&tsirm->S);
    159:   VecDestroy(&tsirm->Alpha);
    160:   VecDestroy(&tsirm->r);
    161:   PetscFree(ksp->data);
    162:   return(0);
    163: }

    165: /*MC
    166:      KSPTSIRM - Implements the two-stage iteration with least-squares residual minimization method.

    169:    Options Database Keys:
    170: +  -ksp_ksp_type <solver> -         the type of the inner solver (GMRES or any of its variants for instance)
    171: .  -ksp_pc_type <preconditioner> - the type of the preconditioner applied to the inner solver
    172: .  -ksp_ksp_max_it <maxits> -      the maximum number of inner iterations (iterations of the inner solver)
    173: .  -ksp_ksp_rtol <tol> -           sets the relative convergence tolerance of the inner solver
    174: .  -ksp_tsirm_cgls <number> -      if 1 use CGLS solver in the minimization step, otherwise use LSQR solver
    175: .  -ksp_tsirm_max_it_ls <maxits> - the maximum number of iterations for the least-squares minimization solver
    176: .  -ksp_tsirm_tol_ls <tol> -       sets the convergence tolerance of the least-squares minimization solver
    177: -  -ksp_tsirm_size_ls <size> -     the number of residuals for the least-squares minimization step

    179:    Level: advanced

    181:    Notes:
    182:     TSIRM is a new two-stage iteration method for solving large sparse linear systems of the form Ax=b. The main idea behind this new
    183:           method is the use a least-squares residual minimization to improve the convergence of Krylov based iterative methods, typically those of GMRES variants.
    184:           The principle of TSIRM algorithm  is to build an outer iteration over a Krylov method, called inner solver, and to frequently store the current residual
    185:           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
    186:           composed by the saved residuals, in order to compute a better solution and to make new iterations if required. The GMRES method , or any of its variants,
    187:           can potentially be used as inner solver. The minimization step consists in solving the least-squares problem min||b-ASa|| to find 'a' which minimizes the
    188:           residuals (b-AS). The minimization step is performed using two solvers of linear least-squares problems: CGLS  or LSQR. A new solution x with
    189:           a minimal residual is computed with x=Sa.

    191:    References:
    192: . 1 R. Couturier, L. Ziane Khodja, and C. Guyeux. TSIRM: A Two-Stage Iteration with least-squares Residual Minimization algorithm to solve large sparse linear systems. In PDSEC 2015, 16th IEEE Int. Workshop on Parallel and Distributed Scientific and Engineering Computing (in conjunction with IPDPS 2015), Hyderabad, India, 2015.

    194:    Contributed by: Lilia Ziane Khodja

    196: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPLGMRES,
    197:            KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
    198:            KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
    199:            KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPSetPCSide()

    201: M*/
    202: PETSC_EXTERN PetscErrorCode KSPCreate_TSIRM(KSP ksp)
    203: {
    205:   KSP_TSIRM      *tsirm;

    208:   PetscNewLog(ksp,&tsirm);
    209:   ksp->data                = (void*)tsirm;
    210:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
    211:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);
    212:   ksp->ops->setup          = KSPSetUp_TSIRM;
    213:   ksp->ops->solve          = KSPSolve_TSIRM;
    214:   ksp->ops->destroy        = KSPDestroy_TSIRM;
    215:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    216:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    217:   ksp->ops->setfromoptions = KSPSetFromOptions_TSIRM;
    218:   ksp->ops->view           = NULL;
    219: #if defined(PETSC_USE_COMPLEX)
    220:   SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"This is not supported for complex numbers");
    221: #else
    222:   return(0);
    223: #endif
    224: }