Actual source code: tsirm.c


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

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

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

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

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

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

 37:   /* Residual and vector Alpha computed in the minimization step */
 38:   PetscCall(MatCreateVecs(tsirm->S, &tsirm->Alpha, &tsirm->r));
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 40: }

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

 57:   PetscFunctionBegin;
 58:   x = ksp->vec_sol; /* Solution vector        */
 59:   b = ksp->vec_rhs; /* Right-hand side vector */

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

 65:   /* Inner solver */
 66:   PetscCall(KSPGetPC(ksp, &pc));
 67:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCKSP, &isksp));
 68:   PetscCheck(isksp, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "PC must be of type PCKSP");
 69:   PetscCall(PCKSPGetKSP(pc, &sub_ksp));
 70:   PetscCall(KSPSetTolerances(sub_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, restart));

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

 87:       /* Build S^T */
 88:       PetscCall(VecGetArray(x, &array));
 89:       PetscCall(MatSetValues(tsirm->S, tsirm->Iend - tsirm->Istart, ind_row, 1, &col, array, INSERT_VALUES));
 90:       PetscCall(VecRestoreArray(x, &array));

 92:       PetscCall(KSPGetResidualNorm(sub_ksp, &norm));
 93:       ksp->rnorm = norm;
 94:       ksp->its++;
 95:       PetscCall(KSPConvergedDefault(ksp, ksp->its, norm, &ksp->reason, ksp->cnvP));
 96:       PetscCall(KSPMonitor(ksp, ksp->its, norm));
 97:     }

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

110:       /* CGLS or LSQR method to minimize the residuals*/
111:       PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_min));
112:       if (tsirm->cgls) {
113:         PetscCall(KSPSetType(ksp_min, KSPCGLS));
114:       } else {
115:         PetscCall(KSPSetType(ksp_min, KSPLSQR));
116:       }
117:       PetscCall(KSPSetOperators(ksp_min, AS, AS));
118:       PetscCall(KSPSetTolerances(ksp_min, tsirm->tol_ls, PETSC_DEFAULT, PETSC_DEFAULT, tsirm->maxiter_ls));
119:       PetscCall(KSPGetPC(ksp_min, &pc_min));
120:       PetscCall(PCSetType(pc_min, PCNONE));
121:       PetscCall(KSPSolve(ksp_min, b, tsirm->Alpha)); /* Find Alpha such that ||AS Alpha = b|| */
122:       PetscCall(KSPDestroy(&ksp_min));
123:       /* Apply minimization */
124:       PetscCall(MatMult(tsirm->S, tsirm->Alpha, x)); /* x = S * Alpha */
125:     }
126:   } while (ksp->its < ksp->max_it && !ksp->reason);
127:   PetscCall(MatDestroy(&AS));
128:   PetscCall(PetscFree(ind_row));
129:   ksp->its = total;
130:   PetscFunctionReturn(PETSC_SUCCESS);
131: }

133: PetscErrorCode KSPSetFromOptions_TSIRM(KSP ksp, PetscOptionItems *PetscOptionsObject)
134: {
135:   KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;

137:   PetscFunctionBegin;
138:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP TSIRM options");
139:   PetscCall(PetscOptionsInt("-ksp_tsirm_cgls", "Method used for the minimization step", "", tsirm->cgls, &tsirm->cgls, NULL)); /*0:LSQR, 1:CGLS*/
140:   PetscCall(PetscOptionsReal("-ksp_tsirm_tol_ls", "Tolerance threshold for the minimization step", "", tsirm->tol_ls, &tsirm->tol_ls, NULL));
141:   PetscCall(PetscOptionsInt("-ksp_tsirm_max_it_ls", "Maximum number of iterations for the minimization step", "", tsirm->maxiter_ls, &tsirm->maxiter_ls, NULL));
142:   PetscCall(PetscOptionsInt("-ksp_tsirm_size_ls", "Number of residuals for minimization", "", tsirm->size_ls, &tsirm->size_ls, NULL));
143:   PetscOptionsHeadEnd();
144:   PetscFunctionReturn(PETSC_SUCCESS);
145: }

147: PetscErrorCode KSPDestroy_TSIRM(KSP ksp)
148: {
149:   KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;

151:   PetscFunctionBegin;
152:   PetscCall(MatDestroy(&tsirm->S));
153:   PetscCall(VecDestroy(&tsirm->Alpha));
154:   PetscCall(VecDestroy(&tsirm->r));
155:   PetscCall(PetscFree(ksp->data));
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

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

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

172:    Level: advanced

174:    Notes:
175:     `KSPTSIRM` is a two-stage iteration method for solving large sparse linear systems of the form Ax=b. The main idea behind this new
176:           method is the use a least-squares residual minimization to improve the convergence of Krylov based iterative methods, typically those of GMRES variants.
177:           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
178:           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
179:           composed by the saved residuals, in order to compute a better solution and to make new iterations if required.
180:           The minimization step consists in solving the least-squares problem min||b-ASa|| to find 'a' which minimizes the
181:           residuals (b-AS). The minimization step is performed using two solvers of linear least-squares problems: `KSPCGLS`  or `KSPLSQR`. A new solution x with
182:           a minimal residual is computed with x=Sa.

184:    References:
185: .  * - 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.
186:    In PDSEC 2015, 16th IEEE Int. Workshop on Parallel and Distributed Scientific and Engineering Computing (in conjunction with IPDPS 2015), Hyderabad, India, 2015.

188:    Contributed by:
189:    Lilia Ziane Khodja

191: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPLGMRES`,
192:           `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
193:           `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
194:           `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPSetPCSide()`
195: M*/
196: PETSC_EXTERN PetscErrorCode KSPCreate_TSIRM(KSP ksp)
197: {
198:   KSP_TSIRM *tsirm;

200:   PetscFunctionBegin;
201:   PetscCall(PetscNew(&tsirm));
202:   ksp->data = (void *)tsirm;
203:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
204:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1));
205:   ksp->ops->setup          = KSPSetUp_TSIRM;
206:   ksp->ops->solve          = KSPSolve_TSIRM;
207:   ksp->ops->destroy        = KSPDestroy_TSIRM;
208:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
209:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
210:   ksp->ops->setfromoptions = KSPSetFromOptions_TSIRM;
211:   ksp->ops->view           = NULL;
212: #if defined(PETSC_USE_COMPLEX)
213:   SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "This is not supported for complex numbers");
214: #else
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: #endif
217: }