Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_ITERATIVE_SOLVER_BASE_H
00026 #define EIGEN_ITERATIVE_SOLVER_BASE_H
00027
00028 namespace Eigen {
00029
00035 template< typename Derived>
00036 class IterativeSolverBase : internal::noncopyable
00037 {
00038 public:
00039 typedef typename internal::traits<Derived>::MatrixType MatrixType;
00040 typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
00041 typedef typename MatrixType::Scalar Scalar;
00042 typedef typename MatrixType::Index Index;
00043 typedef typename MatrixType::RealScalar RealScalar;
00044
00045 public:
00046
00047 Derived& derived() { return *static_cast<Derived*>(this); }
00048 const Derived& derived() const { return *static_cast<const Derived*>(this); }
00049
00051 IterativeSolverBase()
00052 : mp_matrix(0)
00053 {
00054 init();
00055 }
00056
00067 IterativeSolverBase(const MatrixType& A)
00068 {
00069 init();
00070 compute(A);
00071 }
00072
00073 ~IterativeSolverBase() {}
00074
00080 Derived& analyzePattern(const MatrixType& A)
00081 {
00082 m_preconditioner.analyzePattern(A);
00083 m_isInitialized = true;
00084 m_analysisIsOk = true;
00085 m_info = Success;
00086 return derived();
00087 }
00088
00098 Derived& factorize(const MatrixType& A)
00099 {
00100 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00101 mp_matrix = &A;
00102 m_preconditioner.factorize(A);
00103 m_factorizationIsOk = true;
00104 m_info = Success;
00105 return derived();
00106 }
00107
00118 Derived& compute(const MatrixType& A)
00119 {
00120 mp_matrix = &A;
00121 m_preconditioner.compute(A);
00122 m_isInitialized = true;
00123 m_analysisIsOk = true;
00124 m_factorizationIsOk = true;
00125 m_info = Success;
00126 return derived();
00127 }
00128
00130 Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
00132 Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
00133
00135 RealScalar tolerance() const { return m_tolerance; }
00136
00138 Derived& setTolerance(RealScalar tolerance)
00139 {
00140 m_tolerance = tolerance;
00141 return derived();
00142 }
00143
00145 Preconditioner& preconditioner() { return m_preconditioner; }
00146
00148 const Preconditioner& preconditioner() const { return m_preconditioner; }
00149
00151 int maxIterations() const
00152 {
00153 return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations;
00154 }
00155
00157 Derived& setMaxIterations(int maxIters)
00158 {
00159 m_maxIterations = maxIters;
00160 return derived();
00161 }
00162
00164 int iterations() const
00165 {
00166 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
00167 return m_iterations;
00168 }
00169
00171 RealScalar error() const
00172 {
00173 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
00174 return m_error;
00175 }
00176
00181 template<typename Rhs> inline const internal::solve_retval<Derived, Rhs>
00182 solve(const MatrixBase<Rhs>& b) const
00183 {
00184 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00185 eigen_assert(rows()==b.rows()
00186 && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
00187 return internal::solve_retval<Derived, Rhs>(derived(), b.derived());
00188 }
00189
00194 template<typename Rhs>
00195 inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs>
00196 solve(const SparseMatrixBase<Rhs>& b) const
00197 {
00198 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00199 eigen_assert(rows()==b.rows()
00200 && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
00201 return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
00202 }
00203
00205 ComputationInfo info() const
00206 {
00207 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00208 return m_info;
00209 }
00210
00212 template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00213 void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00214 {
00215 eigen_assert(rows()==b.rows());
00216
00217 int rhsCols = b.cols();
00218 int size = b.rows();
00219 Eigen::Matrix<DestScalar,Dynamic,1> tb(size);
00220 Eigen::Matrix<DestScalar,Dynamic,1> tx(size);
00221 for(int k=0; k<rhsCols; ++k)
00222 {
00223 tb = b.col(k);
00224 tx = derived().solve(tb);
00225 dest.col(k) = tx.sparseView(0);
00226 }
00227 }
00228
00229 protected:
00230 void init()
00231 {
00232 m_isInitialized = false;
00233 m_analysisIsOk = false;
00234 m_factorizationIsOk = false;
00235 m_maxIterations = -1;
00236 m_tolerance = NumTraits<Scalar>::epsilon();
00237 }
00238 const MatrixType* mp_matrix;
00239 Preconditioner m_preconditioner;
00240
00241 int m_maxIterations;
00242 RealScalar m_tolerance;
00243
00244 mutable RealScalar m_error;
00245 mutable int m_iterations;
00246 mutable ComputationInfo m_info;
00247 mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk;
00248 };
00249
00250 namespace internal {
00251
00252 template<typename Derived, typename Rhs>
00253 struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs>
00254 : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs>
00255 {
00256 typedef IterativeSolverBase<Derived> Dec;
00257 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00258
00259 template<typename Dest> void evalTo(Dest& dst) const
00260 {
00261 dec().derived()._solve_sparse(rhs(),dst);
00262 }
00263 };
00264
00265 }
00266
00267 }
00268
00269 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H