IterativeSolverBase.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_ITERATIVE_SOLVER_BASE_H
26 #define EIGEN_ITERATIVE_SOLVER_BASE_H
27 
28 namespace Eigen {
29 
35 template< typename Derived>
36 class IterativeSolverBase : internal::noncopyable
37 {
38 public:
39  typedef typename internal::traits<Derived>::MatrixType MatrixType;
40  typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
41  typedef typename MatrixType::Scalar Scalar;
42  typedef typename MatrixType::Index Index;
43  typedef typename MatrixType::RealScalar RealScalar;
44 
45 public:
46 
47  Derived& derived() { return *static_cast<Derived*>(this); }
48  const Derived& derived() const { return *static_cast<const Derived*>(this); }
49 
52  : mp_matrix(0)
53  {
54  init();
55  }
56 
68  {
69  init();
70  compute(A);
71  }
72 
74 
80  Derived& analyzePattern(const MatrixType& A)
81  {
82  m_preconditioner.analyzePattern(A);
83  m_isInitialized = true;
84  m_analysisIsOk = true;
85  m_info = Success;
86  return derived();
87  }
88 
98  Derived& factorize(const MatrixType& A)
99  {
100  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
101  mp_matrix = &A;
102  m_preconditioner.factorize(A);
103  m_factorizationIsOk = true;
104  m_info = Success;
105  return derived();
106  }
107 
118  Derived& compute(const MatrixType& A)
119  {
120  mp_matrix = &A;
121  m_preconditioner.compute(A);
122  m_isInitialized = true;
123  m_analysisIsOk = true;
124  m_factorizationIsOk = true;
125  m_info = Success;
126  return derived();
127  }
128 
130  Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
132  Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
133 
135  RealScalar tolerance() const { return m_tolerance; }
136 
139  {
141  return derived();
142  }
143 
146 
148  const Preconditioner& preconditioner() const { return m_preconditioner; }
149 
151  int maxIterations() const
152  {
153  return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations;
154  }
155 
157  Derived& setMaxIterations(int maxIters)
158  {
159  m_maxIterations = maxIters;
160  return derived();
161  }
162 
164  int iterations() const
165  {
166  eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
167  return m_iterations;
168  }
169 
172  {
173  eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
174  return m_error;
175  }
176 
181  template<typename Rhs> inline const internal::solve_retval<Derived, Rhs>
182  solve(const MatrixBase<Rhs>& b) const
183  {
184  eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
185  eigen_assert(rows()==b.rows()
186  && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
187  return internal::solve_retval<Derived, Rhs>(derived(), b.derived());
188  }
189 
194  template<typename Rhs>
195  inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs>
196  solve(const SparseMatrixBase<Rhs>& b) const
197  {
198  eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
199  eigen_assert(rows()==b.rows()
200  && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
201  return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
202  }
203 
206  {
207  eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
208  return m_info;
209  }
210 
212  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
214  {
215  eigen_assert(rows()==b.rows());
216 
217  int rhsCols = b.cols();
218  int size = b.rows();
221  for(int k=0; k<rhsCols; ++k)
222  {
223  tb = b.col(k);
224  tx = derived().solve(tb);
225  dest.col(k) = tx.sparseView(0);
226  }
227  }
228 
229 protected:
230  void init()
231  {
232  m_isInitialized = false;
233  m_analysisIsOk = false;
234  m_factorizationIsOk = false;
235  m_maxIterations = -1;
237  }
240 
243 
245  mutable int m_iterations;
248 };
249 
250 namespace internal {
251 
252 template<typename Derived, typename Rhs>
253 struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs>
254  : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs>
255 {
256  typedef IterativeSolverBase<Derived> Dec;
258 
259  template<typename Dest> void evalTo(Dest& dst) const
260  {
261  dec().derived()._solve_sparse(rhs(),dst);
262  }
263 };
264 
265 } // end namespace internal
266 
267 } // end namespace Eigen
268 
269 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H