1 #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
2
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;
9
KSPSetUp_TSIRM(KSP ksp)10 static PetscErrorCode KSPSetUp_TSIRM(KSP ksp)
11 {
12 KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;
13
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));
19
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));
25
26 /* Residual and vector Alpha computed in the minimization step */
27 PetscCall(MatCreateVecs(tsirm->S, &tsirm->Alpha, &tsirm->r));
28 PetscFunctionReturn(PETSC_SUCCESS);
29 }
30
KSPSolve_TSIRM(KSP ksp)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;
44
45 PetscFunctionBegin;
46 x = ksp->vec_sol; /* Solution vector */
47 b = ksp->vec_rhs; /* Right-hand side vector */
48
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;
52
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));
58
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;
74
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));
79
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 }
86
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 }
97
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 }
120
KSPSetFromOptions_TSIRM(KSP ksp,PetscOptionItems PetscOptionsObject)121 static PetscErrorCode KSPSetFromOptions_TSIRM(KSP ksp, PetscOptionItems PetscOptionsObject)
122 {
123 KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;
124
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 }
134
KSPDestroy_TSIRM(KSP ksp)135 static PetscErrorCode KSPDestroy_TSIRM(KSP ksp)
136 {
137 KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;
138
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 }
146
147 /*MC
148 KSPTSIRM - Implements the two-stage iteration with least-squares residual minimization method {cite}`couturier2016tsirm`
149
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
159
160 Level: advanced
161
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$.
171
172 Defaults to 30 iterations for the inner solve, use option `-ksp_ksp_max_it <it>` to change it.
173
174 Contributed by:
175 Lilia Ziane Khodja
176
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*/
KSPCreate_TSIRM(KSP ksp)182 PETSC_EXTERN PetscErrorCode KSPCreate_TSIRM(KSP ksp)
183 {
184 KSP_TSIRM *tsirm;
185 PC pc;
186 KSP sub_ksp;
187
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;
208
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 PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "This is not supported for complex numbers");
214 PetscFunctionReturn(PETSC_SUCCESS);
215 }
216