Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #ifndef EIGEN_ITERATIVE_SOLVER_BASE_H
00011 #define EIGEN_ITERATIVE_SOLVER_BASE_H
00012
00013 namespace Eigen {
00014
00020 template< typename Derived>
00021 class IterativeSolverBase : internal::noncopyable
00022 {
00023 public:
00024 typedef typename internal::traits<Derived>::MatrixType MatrixType;
00025 typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
00026 typedef typename MatrixType::Scalar Scalar;
00027 typedef typename MatrixType::Index Index;
00028 typedef typename MatrixType::RealScalar RealScalar;
00029
00030 public:
00031
00032 Derived& derived() { return *static_cast<Derived*>(this); }
00033 const Derived& derived() const { return *static_cast<const Derived*>(this); }
00034
00036 IterativeSolverBase()
00037 : mp_matrix(0)
00038 {
00039 init();
00040 }
00041
00052 IterativeSolverBase(const MatrixType& A)
00053 {
00054 init();
00055 compute(A);
00056 }
00057
00058 ~IterativeSolverBase() {}
00059
00065 Derived& analyzePattern(const MatrixType& A)
00066 {
00067 m_preconditioner.analyzePattern(A);
00068 m_isInitialized = true;
00069 m_analysisIsOk = true;
00070 m_info = Success;
00071 return derived();
00072 }
00073
00083 Derived& factorize(const MatrixType& A)
00084 {
00085 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00086 mp_matrix = &A;
00087 m_preconditioner.factorize(A);
00088 m_factorizationIsOk = true;
00089 m_info = Success;
00090 return derived();
00091 }
00092
00103 Derived& compute(const MatrixType& A)
00104 {
00105 mp_matrix = &A;
00106 m_preconditioner.compute(A);
00107 m_isInitialized = true;
00108 m_analysisIsOk = true;
00109 m_factorizationIsOk = true;
00110 m_info = Success;
00111 return derived();
00112 }
00113
00115 Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
00117 Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
00118
00120 RealScalar tolerance() const { return m_tolerance; }
00121
00123 Derived& setTolerance(const RealScalar& tolerance)
00124 {
00125 m_tolerance = tolerance;
00126 return derived();
00127 }
00128
00130 Preconditioner& preconditioner() { return m_preconditioner; }
00131
00133 const Preconditioner& preconditioner() const { return m_preconditioner; }
00134
00136 int maxIterations() const
00137 {
00138 return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations;
00139 }
00140
00142 Derived& setMaxIterations(int maxIters)
00143 {
00144 m_maxIterations = maxIters;
00145 return derived();
00146 }
00147
00149 int iterations() const
00150 {
00151 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
00152 return m_iterations;
00153 }
00154
00156 RealScalar error() const
00157 {
00158 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
00159 return m_error;
00160 }
00161
00166 template<typename Rhs> inline const internal::solve_retval<Derived, Rhs>
00167 solve(const MatrixBase<Rhs>& b) const
00168 {
00169 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00170 eigen_assert(rows()==b.rows()
00171 && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
00172 return internal::solve_retval<Derived, Rhs>(derived(), b.derived());
00173 }
00174
00179 template<typename Rhs>
00180 inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs>
00181 solve(const SparseMatrixBase<Rhs>& b) const
00182 {
00183 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00184 eigen_assert(rows()==b.rows()
00185 && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
00186 return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
00187 }
00188
00190 ComputationInfo info() const
00191 {
00192 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00193 return m_info;
00194 }
00195
00197 template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00198 void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00199 {
00200 eigen_assert(rows()==b.rows());
00201
00202 int rhsCols = b.cols();
00203 int size = b.rows();
00204 Eigen::Matrix<DestScalar,Dynamic,1> tb(size);
00205 Eigen::Matrix<DestScalar,Dynamic,1> tx(size);
00206 for(int k=0; k<rhsCols; ++k)
00207 {
00208 tb = b.col(k);
00209 tx = derived().solve(tb);
00210 dest.col(k) = tx.sparseView(0);
00211 }
00212 }
00213
00214 protected:
00215 void init()
00216 {
00217 m_isInitialized = false;
00218 m_analysisIsOk = false;
00219 m_factorizationIsOk = false;
00220 m_maxIterations = -1;
00221 m_tolerance = NumTraits<Scalar>::epsilon();
00222 }
00223 const MatrixType* mp_matrix;
00224 Preconditioner m_preconditioner;
00225
00226 int m_maxIterations;
00227 RealScalar m_tolerance;
00228
00229 mutable RealScalar m_error;
00230 mutable int m_iterations;
00231 mutable ComputationInfo m_info;
00232 mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk;
00233 };
00234
00235 namespace internal {
00236
00237 template<typename Derived, typename Rhs>
00238 struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs>
00239 : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs>
00240 {
00241 typedef IterativeSolverBase<Derived> Dec;
00242 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00243
00244 template<typename Dest> void evalTo(Dest& dst) const
00245 {
00246 dec().derived()._solve_sparse(rhs(),dst);
00247 }
00248 };
00249
00250 }
00251
00252 }
00253
00254 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H