Iterative.hxx
1 // Copyright (C) 2003-2009 Marc Duruflé
2 //
3 // This file is part of the linear-algebra library Seldon,
4 // http://seldon.sourceforge.net/.
5 //
6 // Seldon is free software; you can redistribute it and/or modify it under the
7 // terms of the GNU Lesser General Public License as published by the Free
8 // Software Foundation; either version 2.1 of the License, or (at your option)
9 // any later version.
10 //
11 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
12 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
14 // more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public License
17 // along with Seldon. If not, see http://www.gnu.org/licenses/.
18 
19 
20 #ifndef SELDON_FILE_ITERATIVE_HXX
21 
22 namespace Seldon
23 {
25  template<class T>
26  class Preconditioner_Base
27  {
28  public :
29 
30  Preconditioner_Base();
31 
32 #ifdef SELDON_WITH_VIRTUAL
33  virtual ~Preconditioner_Base();
34 
35  virtual void Solve(const VirtualMatrix<T>&, const Vector<T>& r, Vector<T>& z);
36  virtual void TransSolve(const VirtualMatrix<T>&, const Vector<T>& r, Vector<T>& z);
37 
38  virtual void SetInputPreconditioning(const string&, const Vector<string>&);
39  virtual void CopyParameter(const Preconditioner_Base<T>&);
40 #else
41  // solving M z = r
42  template<class Matrix1, class Vector1>
43  void Solve(const Matrix1& A, const Vector1 & r, Vector1 & z);
44 
45  // solving M^t z = r
46  template<class Matrix1, class Vector1>
47  void TransSolve(const Matrix1& A, const Vector1& r, Vector1 & z);
48 #endif
49 
50  };
51 
52 
54 
58  template<class Titer>
59  class Iteration
60  {
61  protected :
62  Titer tolerance;
63  Titer facteur_reste;
64  int max_iter;
65  int nb_iter;
66  int error_code;
67  bool fail_convergence;
68 
74  int print_level;
75  bool init_guess_null;
76  int type_solver;
77  int parameter_restart;
78  int type_preconditioning;
79  string file_name_history;
80 
81  public :
82 
83  Iteration();
84  Iteration(int max_iteration, const Titer& tol);
85  Iteration(const Iteration<Titer>& outer);
86 
87  int GetTypeSolver() const;
88  int GetRestart() const;
89  Titer GetFactor() const;
90  Titer GetTolerance() const;
91  int GetNumberIteration() const;
92 
93  void SetSolver(int type_resolution, int param_restart, int type_prec);
94  void SetRestart(int m);
95  void SetTolerance(Titer stopping_criterion);
96  void SetMaxNumberIteration(int max_iteration);
97  void SetNumberIteration(int nb);
98 
99  void ShowMessages();
100  void ShowFullHistory();
101  void SaveFullHistory(const string& file);
102  void HideMessages();
103 
104  template<class Vector1>
105  int Init(const Vector1& r);
106  bool First() const;
107 
108  bool IsInitGuess_Null() const;
109  void SetInitGuess(bool type);
110 
111  template<class Vector1>
112  bool Finished(const Vector1& r) const;
113  bool Finished(const Titer& r) const;
114 
115  void Fail(int i, const string& s);
116 
117  Iteration& operator ++ (void);
118 
119  int ErrorCode();
120 
121  template<class T1, class Matrix1, class Vector1>
122  void MltAdd(const T1& alpha, const Matrix1& A, const Vector1& x,
123  const T1& beta, Vector1& y);
124 
125  template<class Matrix1, class Vector1>
126  void Mlt(const Matrix1& A, const Vector1& x, Vector1& y);
127 
128  template<class Matrix1, class Vector1>
129  void Mlt(const class_SeldonTrans& trans,
130  const Matrix1& A, const Vector1& x, Vector1& y);
131 
132  };
133 
134  // declarations of all iterative solvers
135 #ifdef SELDON_WITH_VIRTUAL
136  template<class T, class Vector1>
137  int BiCg(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
138  Preconditioner_Base<T>& M,
139  Iteration<typename ClassComplexType<T>::Treal>& iter);
140 
141  template<class T, class Vector1>
142  int BiCgStab(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
143  Preconditioner_Base<T>& M,
144  Iteration<typename ClassComplexType<T>::Treal>& iter);
145 
146  template<class T, class Vector1>
147  int BiCgStabl(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
148  Preconditioner_Base<T>& M,
149  Iteration<typename ClassComplexType<T>::Treal>& iter);
150 
151  template<class T, class Vector1>
152  int BiCgcr(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
153  Preconditioner_Base<T>& M,
154  Iteration<typename ClassComplexType<T>::Treal>& iter);
155 
156  template<class T, class Vector1>
157  int Cg(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
158  Preconditioner_Base<T>& M,
159  Iteration<typename ClassComplexType<T>::Treal>& iter);
160 
161  template<class T, class Vector1>
162  int Cgne(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
163  Preconditioner_Base<T>& M,
164  Iteration<typename ClassComplexType<T>::Treal>& iter);
165 
166  template<class T, class Vector1>
167  int Cgs(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
168  Preconditioner_Base<T>& M,
169  Iteration<typename ClassComplexType<T>::Treal>& iter);
170 
171  template<class T, class Vector1>
172  int CoCg(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
173  Preconditioner_Base<T>& M,
174  Iteration<typename ClassComplexType<T>::Treal>& iter);
175 
176  template<class T, class Vector1>
177  int Gcr(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
178  Preconditioner_Base<T>& M,
179  Iteration<typename ClassComplexType<T>::Treal>& outer);
180 
181  template<class T, class Vector1>
182  int Gmres(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
183  Preconditioner_Base<T>& M,
184  Iteration<typename ClassComplexType<T>::Treal>& outer);
185 
186  template<class T, class Vector1>
187  int Lsqr(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
188  Preconditioner_Base<T>& M,
189  Iteration<typename ClassComplexType<T>::Treal>& iter);
190 
191  template<class T, class Vector1>
192  int MinRes(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
193  Preconditioner_Base<T>& M,
194  Iteration<typename ClassComplexType<T>::Treal>& iter);
195 
196  template<class T, class Vector1>
197  int QCgs(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
198  Preconditioner_Base<T>& M,
199  Iteration<typename ClassComplexType<T>::Treal>& iter);
200 
201  template<class T, class Vector1>
202  int Qmr(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
203  Preconditioner_Base<T>& M,
204  Iteration<typename ClassComplexType<T>::Treal>& iter);
205 
206  template<class T, class Vector1>
207  int QmrSym(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
208  Preconditioner_Base<T>& M,
209  Iteration<typename ClassComplexType<T>::Treal>& iter);
210 
211  template<class T, class Vector1>
212  int Symmlq(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
213  Preconditioner_Base<T>& M,
214  Iteration<typename ClassComplexType<T>::Treal>& iter);
215 
216  template<class T, class Vector1>
217  int TfQmr(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
218  Preconditioner_Base<T>& M,
219  Iteration<typename ClassComplexType<T>::Treal>& iter);
220 
221 #else
222  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
223  int BiCg(const Matrix1& A, Vector1& x, const Vector1& b,
224  Preconditioner& M, Iteration<Titer> & iter);
225 
226  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
227  int BiCgStab(const Matrix1& A, Vector1& x, const Vector1& b,
228  Preconditioner& M, Iteration<Titer> & iter);
229 
230  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
231  int BiCgStabl(const Matrix1& A, Vector1& x, const Vector1& b,
232  Preconditioner& M, Iteration<Titer> & iter);
233 
234  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
235  int BiCgcr(const Matrix1& A, Vector1& x, const Vector1& b,
236  Preconditioner& M, Iteration<Titer> & iter);
237 
238  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
239  int Cg(const Matrix1& A, Vector1& x, const Vector1& b,
240  Preconditioner& M, Iteration<Titer> & iter);
241 
242  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
243  int Cgne(const Matrix1& A, Vector1& x, const Vector1& b,
244  Preconditioner& M, Iteration<Titer> & iter);
245 
246  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
247  int Cgs(const Matrix1& A, Vector1& x, const Vector1& b,
248  Preconditioner& M, Iteration<Titer> & iter);
249 
250  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
251  int CoCg(const Matrix1& A, Vector1& x, const Vector1& b,
252  Preconditioner& M, Iteration<Titer> & iter);
253 
254  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
255  int Gcr(const Matrix1& A, Vector1& x, const Vector1& b,
256  Preconditioner& M, Iteration<Titer> & outer);
257 
258  template <class Titer, class MatrixSparse, class Vector1, class Preconditioner>
259  int Gmres(const MatrixSparse& A, Vector1& x, const Vector1& b,
260  Preconditioner& M, Iteration<Titer> & outer);
261 
262  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
263  int Lsqr(const Matrix1& A, Vector1& x, const Vector1& b,
264  Preconditioner& M, Iteration<Titer> & iter);
265 
266  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
267  int MinRes(const Matrix1& A, Vector1& x, const Vector1& b,
268  Preconditioner& M, Iteration<Titer> & iter);
269 
270  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
271  int QCgs(const Matrix1& A, Vector1& x, const Vector1& b,
272  Preconditioner& M, Iteration<Titer> & iter);
273 
274  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
275  int Qmr(const Matrix1& A, Vector1& x, const Vector1& b,
276  Preconditioner& M, Iteration<Titer> & iter);
277 
278  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
279  int QmrSym(const Matrix1& A, Vector1& x, const Vector1& b,
280  Preconditioner& M, Iteration<Titer> & iter);
281 
282  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
283  int Symmlq(const Matrix1& A, Vector1& x, const Vector1& b,
284  Preconditioner& M, Iteration<Titer> & iter);
285 
286  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
287  int TfQmr(const Matrix1& A, Vector1& x, const Vector1& b,
288  Preconditioner& M, Iteration<Titer> & iter);
289 
290 #endif
291 
292 } // end namespace
293 
294 #define SELDON_FILE_ITERATIVE_HXX
295 #endif
Seldon::Iteration::parameter_restart
int parameter_restart
restart parameter (for Gmres and Gcr)
Definition: Iterative.hxx:77
Seldon::Preconditioner_Base::Preconditioner_Base
Preconditioner_Base()
Default constructor.
Definition: IterativeInline.cxx:27
Seldon::Iteration::First
bool First() const
Returns true if it is the first iteration.
Definition: Iterative.cxx:236
Seldon::Iteration::max_iter
int max_iter
maximum number of iterations
Definition: Iterative.hxx:64
Seldon::Iteration::GetNumberIteration
int GetNumberIteration() const
Returns the number of iterations.
Definition: Iterative.cxx:134
Seldon::Iteration::GetRestart
int GetRestart() const
Returns the restart parameter.
Definition: Iterative.cxx:110
Seldon::QCgs
int QCgs(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves linear system using Quasi-minimized Conjugate Gradient Squared.
Definition: QCgs.cxx:53
Seldon::Iteration::GetTolerance
Titer GetTolerance() const
Returns stopping criterion.
Definition: Iterative.cxx:126
Seldon::TfQmr
int TfQmr(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using Transpose Free Quasi-Minimal Residual.
Definition: TfQmr.cxx:51
Seldon::Iteration::fail_convergence
bool fail_convergence
true if the iterative solver has converged print level
Definition: Iterative.hxx:67
Seldon::Cgs
int Cgs(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves linear system using Conjugate Gradient Squared (CGS)
Definition: Cgs.cxx:52
Seldon::Iteration::type_solver
int type_solver
iterative solver used
Definition: Iterative.hxx:76
Seldon::Vector< T >
Seldon::Iteration::GetTypeSolver
int GetTypeSolver() const
Returns the type of solver.
Definition: Iterative.cxx:102
Seldon::Iteration::SaveFullHistory
void SaveFullHistory(const string &file)
History of residuals is printed in a file.
Definition: Iterative.cxx:201
Seldon::Iteration::SetMaxNumberIteration
void SetMaxNumberIteration(int max_iteration)
Changes the maximum number of iterations.
Definition: Iterative.cxx:169
Seldon::Iteration::Init
int Init(const Vector1 &r)
Initialization with the right hand side.
Definition: Iterative.cxx:218
Seldon::Iteration::operator++
Iteration & operator++(void)
Increment the number of iterations.
Definition: Iterative.cxx:342
Seldon::Iteration::Mlt
void Mlt(const Matrix1 &A, const Vector1 &x, Vector1 &y)
Computes y = A x.
Definition: IterativeInline.cxx:123
Seldon::Cg
int Cg(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using Conjugate Gradient (CG)
Definition: Cg.cxx:52
Seldon::Iteration::ShowMessages
void ShowMessages()
Sets to a normal display (residual each 100 iterations)
Definition: Iterative.cxx:185
Seldon::Iteration::ErrorCode
int ErrorCode()
Returns the error code (if an error occured)
Definition: Iterative.cxx:351
Seldon::BiCgcr
int BiCgcr(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using BiCgCr.
Definition: BiCgcr.cxx:52
Seldon::class_SeldonTrans
Definition: MatrixFlag.hxx:55
Seldon::Iteration::Iteration
Iteration()
Default constructor.
Definition: Iterative.cxx:53
Seldon::Iteration::ShowFullHistory
void ShowFullHistory()
Sets to a complete display (residual each iteration)
Definition: Iterative.cxx:193
Seldon::Gmres
int Gmres(const MatrixSparse &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &outer)
Solves a linear system by using Generalized Minimum Residual (GMRES)
Definition: Gmres.cxx:51
Seldon::Iteration::GetFactor
Titer GetFactor() const
Returns used coefficient to compute relative residual.
Definition: Iterative.cxx:118
Seldon::Iteration::tolerance
Titer tolerance
stopping criterion
Definition: Iterative.hxx:62
Seldon::Iteration::Fail
void Fail(int i, const string &s)
Informs of a failure in the iterative solver.
Definition: Iterative.cxx:330
Seldon::MinRes
int MinRes(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using Minimum Residual (MinRes)
Definition: MinRes.cxx:53
Seldon::Iteration::error_code
int error_code
error code returned by iterative solver
Definition: Iterative.hxx:66
Seldon::Iteration::nb_iter
int nb_iter
number of iterations
Definition: Iterative.hxx:65
Seldon::Iteration::init_guess_null
bool init_guess_null
true if initial guess is null
Definition: Iterative.hxx:75
Seldon::Iteration::MltAdd
void MltAdd(const T1 &alpha, const Matrix1 &A, const Vector1 &x, const T1 &beta, Vector1 &y)
Computes y = beta y + alpha A x.
Definition: IterativeInline.cxx:109
Seldon::Cgne
int Cgne(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system using Conjugate Gradient Normal Equation (CGNE)
Definition: Cgne.cxx:50
Seldon::Symmlq
int Symmlq(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using Symmetric LQ (SymmLQ)
Definition: Symmlq.cxx:52
Seldon::Iteration::SetTolerance
void SetTolerance(Titer stopping_criterion)
Changes the stopping criterion.
Definition: Iterative.cxx:161
Seldon::Preconditioner_Base::TransSolve
void TransSolve(const Matrix1 &A, const Vector1 &r, Vector1 &z)
Solves M^t z = r.
Definition: IterativeInline.cxx:99
Seldon::Iteration::facteur_reste
Titer facteur_reste
inverse of norm of first residual
Definition: Iterative.hxx:63
Seldon::Iteration::SetNumberIteration
void SetNumberIteration(int nb)
Changes the number of iterations.
Definition: Iterative.cxx:177
Seldon::CoCg
int CoCg(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using Conjugate Orthogonal Conjugate Gradient.
Definition: CoCg.cxx:51
Seldon::Qmr
int Qmr(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using Quasi-Minimal Residual (QMR)
Definition: Qmr.cxx:48
Seldon::Lsqr
int Lsqr(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using Least Squares (LSQR)
Definition: Lsqr.cxx:47
Seldon::BiCg
int BiCg(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves a linear system by using BiConjugate Gradient (BICG)
Definition: BiCg.cxx:53
Seldon::Iteration::type_preconditioning
int type_preconditioning
preconditioner used
Definition: Iterative.hxx:78
Seldon::QmrSym
int QmrSym(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Solves linear system using Symmetric Quasi-Minimal Residual (SQMR)
Definition: QmrSym.cxx:48
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::Iteration
Class containing parameters for an iterative resolution.
Definition: Iterative.hxx:59
Seldon::BiCgStabl
int BiCgStabl(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Implements BiConjugate Gradient Stabilized (BICG-STAB(l))
Definition: BiCgStabl.cxx:50
Seldon::Iteration::IsInitGuess_Null
bool IsInitGuess_Null() const
Returns true if the initial guess is null.
Definition: Iterative.cxx:247
Seldon::Gcr
int Gcr(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &outer)
Solves a linear system by using Generalized Conjugate Residual (GCR)
Definition: Gcr.cxx:52
Seldon::Iteration::print_level
int print_level
Definition: Iterative.hxx:74
Seldon::Preconditioner_Base::Solve
void Solve(const Matrix1 &A, const Vector1 &r, Vector1 &z)
Solves M z = r.
Definition: IterativeInline.cxx:87
Seldon::Iteration::SetSolver
void SetSolver(int type_resolution, int param_restart, int type_prec)
Changes the type of solver and preconditioning.
Definition: Iterative.cxx:142
Seldon::Iteration::SetRestart
void SetRestart(int m)
Changes the restart parameter.
Definition: Iterative.cxx:153
Seldon::Iteration::HideMessages
void HideMessages()
Doesn't display any information.
Definition: Iterative.cxx:210
Seldon::VirtualMatrix
Abstract base class for all matrices.
Definition: Matrix_Base.hxx:42
Seldon::Iteration::Finished
bool Finished(const Vector1 &r) const
Returns true if the iterative solver has reached its end.
Definition: Iterative.cxx:264
Seldon::Preconditioner_Base
Base class for preconditioners.
Definition: Iterative.hxx:26
Seldon::BiCgStab
int BiCgStab(const Matrix1 &A, Vector1 &x, const Vector1 &b, Preconditioner &M, Iteration< Titer > &iter)
Implements BiConjugate Gradient Stabilized (BICG-STAB)
Definition: BiCgStab.cxx:50
Seldon::Iteration::SetInitGuess
void SetInitGuess(bool type)
informs if the initial guess is null
Definition: Iterative.cxx:255