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: /* Initialization */
16: #if defined(PETSC_USE_REAL_SINGLE)
17: tsirm->tol_ls = 1e-25;
18: #else
19: tsirm->tol_ls = 1e-50;
20: #endif
21: tsirm->size_ls = 12;
22: tsirm->maxiter_ls = 15;
23: tsirm->cgls = 0;
25: /* Matrix of the system */
26: KSPGetOperators(ksp, &tsirm->A, NULL); /* Matrix of the system */
27: MatGetSize(tsirm->A, &tsirm->size, NULL); /* Size of the system */
28: MatGetOwnershipRange(tsirm->A, &tsirm->Istart, &tsirm->Iend);
30: /* Matrix S of residuals */
31: MatCreate(PETSC_COMM_WORLD, &tsirm->S);
32: MatSetSizes(tsirm->S, tsirm->Iend - tsirm->Istart, PETSC_DECIDE, tsirm->size, tsirm->size_ls);
33: MatSetType(tsirm->S, MATDENSE);
34: MatSetUp(tsirm->S);
36: /* Residual and vector Alpha computed in the minimization step */
37: MatCreateVecs(tsirm->S, &tsirm->Alpha, &tsirm->r);
38: return 0;
39: }
41: PetscErrorCode KSPSolve_TSIRM(KSP ksp)
42: {
43: KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;
44: KSP sub_ksp;
45: PC pc;
46: Mat AS = NULL;
47: Vec x, b;
48: PetscScalar *array;
49: PetscReal norm = 20;
50: PetscInt i, *ind_row, first_iteration = 1, its = 0, total = 0, col = 0;
51: PetscInt restart = 30;
52: KSP ksp_min; /* KSP for minimization */
53: PC pc_min; /* PC for minimization */
54: PetscBool isksp;
56: x = ksp->vec_sol; /* Solution vector */
57: b = ksp->vec_rhs; /* Right-hand side vector */
59: /* Row indexes (these indexes are global) */
60: PetscMalloc1(tsirm->Iend - tsirm->Istart, &ind_row);
61: for (i = 0; i < tsirm->Iend - tsirm->Istart; i++) ind_row[i] = i + tsirm->Istart;
63: /* Inner solver */
64: KSPGetPC(ksp, &pc);
65: PetscObjectTypeCompare((PetscObject)pc, PCKSP, &isksp);
67: PCKSPGetKSP(pc, &sub_ksp);
68: KSPSetTolerances(sub_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, restart);
70: /* previously it seemed good but with SNES it seems not good... */
71: KSP_MatMult(sub_ksp, tsirm->A, x, tsirm->r);
72: VecAXPY(tsirm->r, -1, b);
73: VecNorm(tsirm->r, NORM_2, &norm);
74: KSPCheckNorm(ksp, norm);
75: ksp->its = 0;
76: KSPConvergedDefault(ksp, ksp->its, norm, &ksp->reason, ksp->cnvP);
77: KSPSetInitialGuessNonzero(sub_ksp, PETSC_TRUE);
78: do {
79: for (col = 0; col < tsirm->size_ls && ksp->reason == 0; col++) {
80: /* Solve (inner iteration) */
81: KSPSolve(sub_ksp, b, x);
82: KSPGetIterationNumber(sub_ksp, &its);
83: total += its;
85: /* Build S^T */
86: VecGetArray(x, &array);
87: MatSetValues(tsirm->S, tsirm->Iend - tsirm->Istart, ind_row, 1, &col, array, INSERT_VALUES);
88: VecRestoreArray(x, &array);
90: KSPGetResidualNorm(sub_ksp, &norm);
91: ksp->rnorm = norm;
92: ksp->its++;
93: KSPConvergedDefault(ksp, ksp->its, norm, &ksp->reason, ksp->cnvP);
94: KSPMonitor(ksp, ksp->its, norm);
95: }
97: /* Minimization step */
98: if (!ksp->reason) {
99: MatAssemblyBegin(tsirm->S, MAT_FINAL_ASSEMBLY);
100: MatAssemblyEnd(tsirm->S, MAT_FINAL_ASSEMBLY);
101: if (first_iteration) {
102: MatMatMult(tsirm->A, tsirm->S, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &AS);
103: first_iteration = 0;
104: } else {
105: MatMatMult(tsirm->A, tsirm->S, MAT_REUSE_MATRIX, PETSC_DEFAULT, &AS);
106: }
108: /* CGLS or LSQR method to minimize the residuals*/
109: KSPCreate(PETSC_COMM_WORLD, &ksp_min);
110: if (tsirm->cgls) {
111: KSPSetType(ksp_min, KSPCGLS);
112: } else {
113: KSPSetType(ksp_min, KSPLSQR);
114: }
115: KSPSetOperators(ksp_min, AS, AS);
116: KSPSetTolerances(ksp_min, tsirm->tol_ls, PETSC_DEFAULT, PETSC_DEFAULT, tsirm->maxiter_ls);
117: KSPGetPC(ksp_min, &pc_min);
118: PCSetType(pc_min, PCNONE);
119: KSPSolve(ksp_min, b, tsirm->Alpha); /* Find Alpha such that ||AS Alpha = b|| */
120: KSPDestroy(&ksp_min);
121: /* Apply minimization */
122: MatMult(tsirm->S, tsirm->Alpha, x); /* x = S * Alpha */
123: }
124: } while (ksp->its < ksp->max_it && !ksp->reason);
125: MatDestroy(&AS);
126: PetscFree(ind_row);
127: ksp->its = total;
128: return 0;
129: }
131: PetscErrorCode KSPSetFromOptions_TSIRM(KSP ksp, PetscOptionItems *PetscOptionsObject)
132: {
133: KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;
135: PetscOptionsHeadBegin(PetscOptionsObject, "KSP TSIRM options");
136: PetscOptionsInt("-ksp_tsirm_cgls", "Method used for the minimization step", "", tsirm->cgls, &tsirm->cgls, NULL); /*0:LSQR, 1:CGLS*/
137: PetscOptionsReal("-ksp_tsirm_tol_ls", "Tolerance threshold for the minimization step", "", tsirm->tol_ls, &tsirm->tol_ls, NULL);
138: PetscOptionsInt("-ksp_tsirm_max_it_ls", "Maximum number of iterations for the minimization step", "", tsirm->maxiter_ls, &tsirm->maxiter_ls, NULL);
139: PetscOptionsInt("-ksp_tsirm_size_ls", "Number of residuals for minimization", "", tsirm->size_ls, &tsirm->size_ls, NULL);
140: PetscOptionsHeadEnd();
141: return 0;
142: }
144: PetscErrorCode KSPDestroy_TSIRM(KSP ksp)
145: {
146: KSP_TSIRM *tsirm = (KSP_TSIRM *)ksp->data;
148: MatDestroy(&tsirm->S);
149: VecDestroy(&tsirm->Alpha);
150: VecDestroy(&tsirm->r);
151: PetscFree(ksp->data);
152: return 0;
153: }
155: /*MC
156: KSPTSIRM - Implements the two-stage iteration with least-squares residual minimization method.
158: Options Database Keys:
159: + -ksp_ksp_type <solver> - the type of the inner solver (GMRES or any of its variants for instance)
160: . -ksp_pc_type <preconditioner> - the type of the preconditioner applied to the inner solver
161: . -ksp_ksp_max_it <maxits> - the maximum number of inner iterations (iterations of the inner solver)
162: . -ksp_ksp_rtol <tol> - sets the relative convergence tolerance of the inner solver
163: . -ksp_tsirm_cgls <number> - if 1 use CGLS solver in the minimization step, otherwise use LSQR solver
164: . -ksp_tsirm_max_it_ls <maxits> - the maximum number of iterations for the least-squares minimization solver
165: . -ksp_tsirm_tol_ls <tol> - sets the convergence tolerance of the least-squares minimization solver
166: - -ksp_tsirm_size_ls <size> - the number of residuals for the least-squares minimization step
168: Level: advanced
170: Notes:
171: `KSPTSIRM` is a two-stage iteration method for solving large sparse linear systems of the form Ax=b. The main idea behind this new
172: method is the use a least-squares residual minimization to improve the convergence of Krylov based iterative methods, typically those of GMRES variants.
173: 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
174: 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
175: composed by the saved residuals, in order to compute a better solution and to make new iterations if required.
176: The minimization step consists in solving the least-squares problem min||b-ASa|| to find 'a' which minimizes the
177: residuals (b-AS). The minimization step is performed using two solvers of linear least-squares problems: `KSPCGLS` or `KSPLSQR`. A new solution x with
178: a minimal residual is computed with x=Sa.
180: References:
181: . * - 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.
182: In PDSEC 2015, 16th IEEE Int. Workshop on Parallel and Distributed Scientific and Engineering Computing (in conjunction with IPDPS 2015), Hyderabad, India, 2015.
184: Contributed by:
185: Lilia Ziane Khodja
187: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPLGMRES`,
188: `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
189: `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
190: `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPSetPCSide()`
191: M*/
192: PETSC_EXTERN PetscErrorCode KSPCreate_TSIRM(KSP ksp)
193: {
194: KSP_TSIRM *tsirm;
196: PetscNew(&tsirm);
197: ksp->data = (void *)tsirm;
198: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
199: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1);
200: ksp->ops->setup = KSPSetUp_TSIRM;
201: ksp->ops->solve = KSPSolve_TSIRM;
202: ksp->ops->destroy = KSPDestroy_TSIRM;
203: ksp->ops->buildsolution = KSPBuildSolutionDefault;
204: ksp->ops->buildresidual = KSPBuildResidualDefault;
205: ksp->ops->setfromoptions = KSPSetFromOptions_TSIRM;
206: ksp->ops->view = NULL;
207: #if defined(PETSC_USE_COMPLEX)
208: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "This is not supported for complex numbers");
209: #else
210: return 0;
211: #endif
212: }