GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-01 08:14:10 +08:00
Code Issues Packages Projects Releases Wiki Activity

Make unsupport sparse solvers use SparseSolverBase

Browse Source
This commit is contained in:
Gael Guennebaud 2014-09-01 17:21:47 +02:00
parent b3d63b4db2
commit b051bbd64f
5 changed files with 55 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
unsupported/Eigen/src/IterativeSolvers/DGMRES.h
View File

@ -108,6 +108,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
using Base::m_isInitialized; using Base::m_isInitialized;
using Base::m_tolerance; using Base::m_tolerance;
public: public:
using Base::_solve_impl;
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
@ -138,6 +139,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
~DGMRES() {} ~DGMRES() {}
#ifndef EIGEN_TEST_EVALUATORS
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution. * \a x0 as an initial solution.
* *
@ -153,10 +155,11 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
return internal::solve_retval_with_guess return internal::solve_retval_with_guess
<DGMRES, Rhs, Guess>(*this, b.derived(), x0); <DGMRES, Rhs, Guess>(*this, b.derived(), x0);
} }
#endif
/** \internal */ /** \internal */
template<typename Rhs,typename Dest> template<typename Rhs,typename Dest>
void _solveWithGuess(const Rhs& b, Dest& x) const void _solve_with_guess_impl(const Rhs& b, Dest& x) const
{ {
bool failed = false; bool failed = false;
for(int j=0; j<b.cols(); ++j) for(int j=0; j<b.cols(); ++j)
@ -175,10 +178,10 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
/** \internal */ /** \internal */
template<typename Rhs,typename Dest> template<typename Rhs,typename Dest>
void _solve(const Rhs& b, Dest& x) const void _solve_impl(const Rhs& b, MatrixBase<Dest>& x) const
{ {
x = b; x = b;
_solveWithGuess(b,x); _solve_with_guess_impl(b,x.derived());
} }
/** /**
* Get the restart value * Get the restart value
@ -522,6 +525,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x,
return 0; return 0;
} }
#ifndef EIGEN_TEST_EVALUATORS
namespace internal { namespace internal {
template<typename _MatrixType, typename _Preconditioner, typename Rhs> template<typename _MatrixType, typename _Preconditioner, typename Rhs>
@ -533,10 +537,11 @@ struct solve_retval<DGMRES<_MatrixType, _Preconditioner>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const template<typename Dest> void evalTo(Dest& dst) const
{ {
dec()._solve(rhs(),dst); dec()._solve_impl(rhs(),dst);
} }
}; };
} // end namespace internal } // end namespace internal
#endif
} // end namespace Eigen } // end namespace Eigen
#endif #endif

15
unsupported/Eigen/src/IterativeSolvers/GMRES.h
View File

@ -281,6 +281,7 @@ private:
int m_restart; int m_restart;
public: public:
using Base::_solve_impl;
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
@ -315,6 +316,7 @@ public:
*/ */
void set_restart(const int restart) { m_restart=restart; } void set_restart(const int restart) { m_restart=restart; }
#ifndef EIGEN_TEST_EVALUATORS
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution. * \a x0 as an initial solution.
* *
@ -330,10 +332,11 @@ public:
return internal::solve_retval_with_guess return internal::solve_retval_with_guess
<GMRES, Rhs, Guess>(*this, b.derived(), x0); <GMRES, Rhs, Guess>(*this, b.derived(), x0);
} }
#endif
/** \internal */ /** \internal */
template<typename Rhs,typename Dest> template<typename Rhs,typename Dest>
void _solveWithGuess(const Rhs& b, Dest& x) const void _solve_with_guess_impl(const Rhs& b, Dest& x) const
{ {
bool failed = false; bool failed = false;
for(int j=0; j<b.cols(); ++j) for(int j=0; j<b.cols(); ++j)
@ -353,18 +356,18 @@ public:
/** \internal */ /** \internal */
template<typename Rhs,typename Dest> template<typename Rhs,typename Dest>
void _solve(const Rhs& b, Dest& x) const void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const
{ {
x = b; x = b;
if(x.squaredNorm() == 0) return; // Check Zero right hand side if(x.squaredNorm() == 0) return; // Check Zero right hand side
_solveWithGuess(b,x); _solve_with_guess_impl(b,x.derived());
} }
protected: protected:
}; };
#ifndef EIGEN_TEST_EVALUATORS
namespace internal { namespace internal {
template<typename _MatrixType, typename _Preconditioner, typename Rhs> template<typename _MatrixType, typename _Preconditioner, typename Rhs>
@ -376,12 +379,12 @@ struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const template<typename Dest> void evalTo(Dest& dst) const
{ {
dec()._solve(rhs(),dst); dec()._solve_impl(rhs(),dst);
} }
}; };
} // end namespace internal } // end namespace internal
#endif
} // end namespace Eigen } // end namespace Eigen
#endif // EIGEN_GMRES_H #endif // EIGEN_GMRES_H

15
unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
View File

@ -27,8 +27,11 @@ namespace Eigen {
*/ */
template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> > template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> >
class IncompleteCholesky : internal::noncopyable class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
{ {
protected:
typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
using Base::m_isInitialized;
public: public:
typedef SparseMatrix<Scalar,ColMajor> MatrixType; typedef SparseMatrix<Scalar,ColMajor> MatrixType;
typedef _OrderingType OrderingType; typedef _OrderingType OrderingType;
@ -89,7 +92,7 @@ class IncompleteCholesky : internal::noncopyable
} }
template<typename Rhs, typename Dest> template<typename Rhs, typename Dest>
void _solve(const Rhs& b, Dest& x) const void _solve_impl(const Rhs& b, Dest& x) const
{ {
eigen_assert(m_factorizationIsOk && "factorize() should be called first"); eigen_assert(m_factorizationIsOk && "factorize() should be called first");
if (m_perm.rows() == b.rows()) if (m_perm.rows() == b.rows())
@ -103,6 +106,8 @@ class IncompleteCholesky : internal::noncopyable
x = m_perm * x; x = m_perm * x;
x = m_scal.asDiagonal() * x; x = m_scal.asDiagonal() * x;
} }
#ifndef EIGEN_TEST_EVALUATORS
template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs> template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs>
solve(const MatrixBase<Rhs>& b) const solve(const MatrixBase<Rhs>& b) const
{ {
@ -112,13 +117,14 @@ class IncompleteCholesky : internal::noncopyable
&& "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b"); && "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b");
return internal::solve_retval<IncompleteCholesky, Rhs>(*this, b.derived()); return internal::solve_retval<IncompleteCholesky, Rhs>(*this, b.derived());
} }
#endif
protected: protected:
SparseMatrix<Scalar,ColMajor> m_L; // The lower part stored in CSC SparseMatrix<Scalar,ColMajor> m_L; // The lower part stored in CSC
ScalarType m_scal; // The vector for scaling the matrix ScalarType m_scal; // The vector for scaling the matrix
Scalar m_shift; //The initial shift parameter Scalar m_shift; //The initial shift parameter
bool m_analysisIsOk; bool m_analysisIsOk;
bool m_factorizationIsOk; bool m_factorizationIsOk;
bool m_isInitialized;
ComputationInfo m_info; ComputationInfo m_info;
PermutationType m_perm; PermutationType m_perm;
@ -256,6 +262,8 @@ inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const Idx
listCol[rowIdx(jk)].push_back(col); listCol[rowIdx(jk)].push_back(col);
} }
} }
#ifndef EIGEN_TEST_EVALUATORS
namespace internal { namespace internal {
template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs> template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs>
@ -272,6 +280,7 @@ struct solve_retval<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs>
}; };
} // end namespace internal } // end namespace internal
#endif
} // end namespace Eigen } // end namespace Eigen

19
unsupported/Eigen/src/IterativeSolvers/IncompleteLU.h
View File

@ -13,8 +13,12 @@
namespace Eigen { namespace Eigen {
template <typename _Scalar> template <typename _Scalar>
class IncompleteLU class IncompleteLU : public SparseSolverBase<IncompleteLU<_Scalar> >
{ {
protected:
typedef SparseSolverBase<IncompleteLU<_Scalar> > Base;
using Base::m_isInitialized;
typedef _Scalar Scalar; typedef _Scalar Scalar;
typedef Matrix<Scalar,Dynamic,1> Vector; typedef Matrix<Scalar,Dynamic,1> Vector;
typedef typename Vector::Index Index; typedef typename Vector::Index Index;
@ -23,10 +27,10 @@ class IncompleteLU
public: public:
typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
IncompleteLU() : m_isInitialized(false) {} IncompleteLU() {}
template<typename MatrixType> template<typename MatrixType>
IncompleteLU(const MatrixType& mat) : m_isInitialized(false) IncompleteLU(const MatrixType& mat)
{ {
compute(mat); compute(mat);
} }
@ -71,12 +75,13 @@ class IncompleteLU
} }
template<typename Rhs, typename Dest> template<typename Rhs, typename Dest>
void _solve(const Rhs& b, Dest& x) const void _solve_impl(const Rhs& b, Dest& x) const
{ {
x = m_lu.template triangularView<UnitLower>().solve(b); x = m_lu.template triangularView<UnitLower>().solve(b);
x = m_lu.template triangularView<Upper>().solve(x); x = m_lu.template triangularView<Upper>().solve(x);
} }
#ifndef EIGEN_TEST_EVALUATORS
template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs> template<typename Rhs> inline const internal::solve_retval<IncompleteLU, Rhs>
solve(const MatrixBase<Rhs>& b) const solve(const MatrixBase<Rhs>& b) const
{ {
@ -85,12 +90,13 @@ class IncompleteLU
&& "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b"); && "IncompleteLU::solve(): invalid number of rows of the right hand side matrix b");
return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived()); return internal::solve_retval<IncompleteLU, Rhs>(*this, b.derived());
} }
#endif
protected: protected:
FactorType m_lu; FactorType m_lu;
bool m_isInitialized;
}; };
#ifndef EIGEN_TEST_EVALUATORS
namespace internal { namespace internal {
template<typename _MatrixType, typename Rhs> template<typename _MatrixType, typename Rhs>
@ -102,11 +108,12 @@ struct solve_retval<IncompleteLU<_MatrixType>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const template<typename Dest> void evalTo(Dest& dst) const
{ {
dec()._solve(rhs(),dst); dec()._solve_impl(rhs(),dst);
} }
}; };
} // end namespace internal } // end namespace internal
#endif
} // end namespace Eigen } // end namespace Eigen

15
unsupported/Eigen/src/IterativeSolvers/MINRES.h
View File

@ -2,7 +2,7 @@
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2012 Giacomo Po <gpo@ucla.edu> // Copyright (C) 2012 Giacomo Po <gpo@ucla.edu>
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
// //
// This Source Code Form is subject to the terms of the Mozilla // This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed // Public License v. 2.0. If a copy of the MPL was not distributed
@ -217,6 +217,7 @@ namespace Eigen {
using Base::m_info; using Base::m_info;
using Base::m_isInitialized; using Base::m_isInitialized;
public: public:
using Base::_solve_impl;
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
@ -245,6 +246,7 @@ namespace Eigen {
/** Destructor. */ /** Destructor. */
~MINRES(){} ~MINRES(){}
#ifndef EIGEN_TEST_EVALUATORS
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution. * \a x0 as an initial solution.
* *
@ -260,10 +262,11 @@ namespace Eigen {
return internal::solve_retval_with_guess return internal::solve_retval_with_guess
<MINRES, Rhs, Guess>(*this, b.derived(), x0); <MINRES, Rhs, Guess>(*this, b.derived(), x0);
} }
#endif
/** \internal */ /** \internal */
template<typename Rhs,typename Dest> template<typename Rhs,typename Dest>
void _solveWithGuess(const Rhs& b, Dest& x) const void _solve_with_guess_impl(const Rhs& b, Dest& x) const
{ {
m_iterations = Base::maxIterations(); m_iterations = Base::maxIterations();
m_error = Base::m_tolerance; m_error = Base::m_tolerance;
@ -284,16 +287,17 @@ namespace Eigen {
/** \internal */ /** \internal */
template<typename Rhs,typename Dest> template<typename Rhs,typename Dest>
void _solve(const Rhs& b, Dest& x) const void _solve_impl(const Rhs& b, MatrixBase<Dest> &x) const
{ {
x.setZero(); x.setZero();
_solveWithGuess(b,x); _solve_with_guess_impl(b,x.derived());
} }
protected: protected:
}; };
#ifndef EIGEN_TEST_EVALUATORS
namespace internal { namespace internal {
template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs> template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
@ -305,11 +309,12 @@ namespace Eigen {
template<typename Dest> void evalTo(Dest& dst) const template<typename Dest> void evalTo(Dest& dst) const
{ {
dec()._solve(rhs(),dst); dec()._solve_impl(rhs(),dst);
} }
}; };
} // end namespace internal } // end namespace internal
#endif
} // end namespace Eigen } // end namespace Eigen