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: }