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*/
110: KSPCreate(PETSC_COMM_WORLD,&ksp_min);
111: if (tsirm->cgls) {
112: KSPSetType(ksp_min,KSPCGLS);
113: } else {
114: KSPSetType(ksp_min,KSPLSQR);
115: }
116: KSPSetOperators(ksp_min,AS,AS);
117: KSPSetTolerances(ksp_min,tsirm->tol_ls,PETSC_DEFAULT,PETSC_DEFAULT,tsirm->maxiter_ls);
118: KSPGetPC(ksp_min,&pc_min);
119: PCSetType(pc_min,PCNONE);
120: KSPSolve(ksp_min,b,tsirm->Alpha); /* Find Alpha such that ||AS Alpha = b|| */
121: KSPDestroy(&ksp_min);
122: /* Apply minimization */
123: MatMult(tsirm->S,tsirm->Alpha,x); /* x = S * Alpha */
124: }
125: } while (ksp->its<ksp->max_it && !ksp->reason);
126: MatDestroy(&AS);
127: PetscFree(ind_row);
128: ksp->its = total;
129: return 0;
130: }
132: PetscErrorCode KSPSetFromOptions_TSIRM(PetscOptionItems *PetscOptionsObject,KSP ksp)
133: {
134: KSP_TSIRM *tsirm = (KSP_TSIRM*)ksp->data;
136: PetscOptionsHead(PetscOptionsObject,"KSP TSIRM options");
137: PetscOptionsInt("-ksp_tsirm_cgls","Method used for the minimization step","",tsirm->cgls,&tsirm->cgls,NULL); /*0:LSQR, 1:CGLS*/
138: PetscOptionsReal("-ksp_tsirm_tol_ls","Tolerance threshold for the minimization step","",tsirm->tol_ls,&tsirm->tol_ls,NULL);
139: PetscOptionsInt("-ksp_tsirm_max_it_ls","Maximum number of iterations for the minimization step","",tsirm->maxiter_ls,&tsirm->maxiter_ls,NULL);
140: PetscOptionsInt("-ksp_tsirm_size_ls","Number of residuals for minimization","",tsirm->size_ls,&tsirm->size_ls,NULL);
141: PetscOptionsTail();
142: return 0;
143: }
145: PetscErrorCode KSPDestroy_TSIRM(KSP ksp)
146: {
147: KSP_TSIRM *tsirm = (KSP_TSIRM*)ksp->data;
149: MatDestroy(&tsirm->S);
150: VecDestroy(&tsirm->Alpha);
151: VecDestroy(&tsirm->r);
152: PetscFree(ksp->data);
153: return 0;
154: }
156: /*MC
157: KSPTSIRM - Implements the two-stage iteration with least-squares residual minimization method.
159: Options Database Keys:
160: + -ksp_ksp_type <solver> - the type of the inner solver (GMRES or any of its variants for instance)
161: . -ksp_pc_type <preconditioner> - the type of the preconditioner applied to the inner solver
162: . -ksp_ksp_max_it <maxits> - the maximum number of inner iterations (iterations of the inner solver)
163: . -ksp_ksp_rtol <tol> - sets the relative convergence tolerance of the inner solver
164: . -ksp_tsirm_cgls <number> - if 1 use CGLS solver in the minimization step, otherwise use LSQR solver
165: . -ksp_tsirm_max_it_ls <maxits> - the maximum number of iterations for the least-squares minimization solver
166: . -ksp_tsirm_tol_ls <tol> - sets the convergence tolerance of the least-squares minimization solver
167: - -ksp_tsirm_size_ls <size> - the number of residuals for the least-squares minimization step
169: Level: advanced
171: Notes:
172: 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
173: method is the use a least-squares residual minimization to improve the convergence of Krylov based iterative methods, typically those of GMRES variants.
174: 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
175: 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
176: 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,
177: 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
178: residuals (b-AS). The minimization step is performed using two solvers of linear least-squares problems: CGLS or LSQR. A new solution x with
179: a minimal residual is computed with x=Sa.
181: References:
182: . * - 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.
184: Contributed by: Lilia Ziane Khodja
186: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPLGMRES,
187: KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
188: KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
189: KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov(), KSPSetPCSide()
191: M*/
192: PETSC_EXTERN PetscErrorCode KSPCreate_TSIRM(KSP ksp)
193: {
194: KSP_TSIRM *tsirm;
196: PetscNewLog(ksp,&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: }