devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Refactoring of sparse solvers through a SparseSolverBase class and usage of the Solve<> expression. Introduce a SolveWithGuess expression on top of Solve.

Browse Source
This commit is contained in:
Gael Guennebaud
2014-09-01 15:00:19 +02:00
parent bc065c75d2
commit 85c7659574
20 changed files with 565 additions and 131 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h
View File

@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -80,12 +80,14 @@ class DiagonalPreconditioner
return factorize(mat);
}
/** \internal */
template<typename Rhs, typename Dest>
void _solve(const Rhs& b, Dest& x) const
void _solve_impl(const Rhs& b, Dest& x) const
{
x = m_invdiag.array() * b.array() ;
}
#ifndef EIGEN_TEST_EVALUATORS
template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
@@ -94,12 +96,23 @@ class DiagonalPreconditioner
&& "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived());
}
#else
template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
eigen_assert(m_invdiag.size()==b.rows()
&& "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived());
}
#endif
protected:
Vector m_invdiag;
bool m_isInitialized;
};
#ifndef EIGEN_TEST_EVALUATORS
namespace internal {
template<typename _MatrixType, typename Rhs>
@@ -111,11 +124,12 @@ struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dec()._solve(rhs(),dst);
dec()._solve_impl(rhs(),dst);
}
};
}
#endif // EIGEN_TEST_EVALUATORS
/** \ingroup IterativeLinearSolvers_Module
* \brief A naive preconditioner which approximates any matrix as the identity matrix

41
Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
View File

@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
@@ -181,7 +181,8 @@ public:
BiCGSTAB(const MatrixType& A) : Base(A) {}
~BiCGSTAB() {}
#ifndef EIGEN_TEST_EVALUATORS
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution.
*
@@ -197,10 +198,26 @@ public:
return internal::solve_retval_with_guess
<BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0);
}
#else
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution.
*
* \sa compute()
*/
template<typename Rhs,typename Guess>
inline const SolveWithGuess<BiCGSTAB, Rhs, Guess>
solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
{
eigen_assert(m_isInitialized && "BiCGSTAB is not initialized.");
eigen_assert(Base::rows()==b.rows()
&& "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b");
return SolveWithGuess<BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0);
}
#endif
/** \internal */
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;
for(int j=0; j<b.cols(); ++j)
@@ -219,22 +236,23 @@ public:
}
/** \internal */
using Base::_solve_impl;
template<typename Rhs,typename Dest>
void _solve(const Rhs& b, Dest& x) const
void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const
{
// x.setZero();
x = b;
_solveWithGuess(b,x);
// x.setZero();
x = b;
_solve_with_guess_impl(b,x);
}
protected:
};
#ifndef EIGEN_TEST_EVALUATORS
namespace internal {
template<typename _MatrixType, typename _Preconditioner, typename Rhs>
template<typename _MatrixType, typename _Preconditioner, typename Rhs>
struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
: solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
{
@@ -243,11 +261,12 @@ struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dec()._solve(rhs(),dst);
dec()._solve_impl(rhs(),dst);
}
};
} // end namespace internal
#endif
} // end namespace Eigen

31
Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
View File

@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -193,6 +193,7 @@ public:
~ConjugateGradient() {}
#ifndef EIGEN_TEST_EVALUATORS
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution.
*
@@ -208,10 +209,26 @@ public:
return internal::solve_retval_with_guess
<ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0);
}
#else
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution.
*
* \sa compute()
*/
template<typename Rhs,typename Guess>
inline const SolveWithGuess<ConjugateGradient, Rhs, Guess>
solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
{
eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
eigen_assert(Base::rows()==b.rows()
&& "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
return SolveWithGuess<ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0);
}
#endif
/** \internal */
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_error = Base::m_tolerance;
@@ -231,18 +248,19 @@ public:
}
/** \internal */
using Base::_solve_impl;
template<typename Rhs,typename Dest>
void _solve(const Rhs& b, Dest& x) const
void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const
{
x.setOnes();
_solveWithGuess(b,x);
_solve_with_guess_impl(b.derived(),x);
}
protected:
};
#ifndef EIGEN_TEST_EVALUATORS
namespace internal {
template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
@@ -254,11 +272,12 @@ struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
dec()._solve(rhs(),dst);
dec()._solve_impl(rhs(),dst);
}
};
} // end namespace internal
#endif
} // end namespace Eigen

20
Eigen/src/IterativeLinearSolvers/IncompleteLUT.h
View File

@@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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
@@ -158,7 +159,11 @@ class IncompleteLUT : internal::noncopyable
void setFillfactor(int fillfactor);
template<typename Rhs, typename Dest>
#ifndef EIGEN_TEST_EVALUATORS
void _solve(const Rhs& b, Dest& x) const
#else
void _solve_impl(const Rhs& b, Dest& x) const
#endif
{
x = m_Pinv * b;
x = m_lu.template triangularView<UnitLower>().solve(x);
@@ -166,14 +171,25 @@ class IncompleteLUT : internal::noncopyable
x = m_P * x;
}
#ifndef EIGEN_TEST_EVALUATORS
template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs>
solve(const MatrixBase<Rhs>& b) const
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
eigen_assert(cols()==b.rows()
&& "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived());
}
#else
template<typename Rhs> inline const Solve<IncompleteLUT, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
eigen_assert(cols()==b.rows()
&& "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
return Solve<IncompleteLUT, Rhs>(*this, b.derived());
}
#endif
protected:
@@ -445,6 +461,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
m_info = Success;
}
#ifndef EIGEN_TEST_EVALUATORS
namespace internal {
template<typename _MatrixType, typename Rhs>
@@ -461,6 +478,7 @@ struct solve_retval<IncompleteLUT<_MatrixType>, Rhs>
};
} // end namespace internal
#endif // EIGEN_TEST_EVALUATORS
} // end namespace Eigen

39
Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h
View File

@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -18,8 +18,12 @@ namespace Eigen {
* \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
*/
template< typename Derived>
class IterativeSolverBase : internal::noncopyable
class IterativeSolverBase : public SparseSolverBase<Derived>
{
protected:
typedef SparseSolverBase<Derived> Base;
using Base::m_isInitialized;
public:
typedef typename internal::traits<Derived>::MatrixType MatrixType;
typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
@@ -29,8 +33,7 @@ public:
public:
Derived& derived() { return *static_cast<Derived*>(this); }
const Derived& derived() const { return *static_cast<const Derived*>(this); }
using Base::derived;
/** Default constructor. */
IterativeSolverBase()
@@ -159,6 +162,7 @@ public:
return m_error;
}
#ifndef EIGEN_TEST_EVALUATORS
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
*
* \sa compute()
@@ -171,7 +175,9 @@ public:
&& "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
return internal::solve_retval<Derived, Rhs>(derived(), b.derived());
}
#endif
#ifndef EIGEN_TEST_EVALUATORS
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
*
* \sa compute()
@@ -185,6 +191,7 @@ public:
&& "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
}
#endif // EIGEN_TEST_EVALUATORS
/** \returns Success if the iterations converged, and NoConvergence otherwise. */
ComputationInfo info() const
@@ -193,6 +200,7 @@ public:
return m_info;
}
#ifndef EIGEN_TEST_EVALUATORS
/** \internal */
template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
@@ -210,6 +218,25 @@ public:
dest.col(k) = tx.sparseView(0);
}
}
#else
/** \internal */
template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
void _solve_impl(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
{
eigen_assert(rows()==b.rows());
int rhsCols = b.cols();
int size = b.rows();
Eigen::Matrix<DestScalar,Dynamic,1> tb(size);
Eigen::Matrix<DestScalar,Dynamic,1> tx(size);
for(int k=0; k<rhsCols; ++k)
{
tb = b.col(k);
tx = derived().solve(tb);
dest.col(k) = tx.sparseView(0);
}
}
#endif // EIGEN_TEST_EVALUATORS
protected:
void init()
@@ -229,9 +256,10 @@ protected:
mutable RealScalar m_error;
mutable int m_iterations;
mutable ComputationInfo m_info;
mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk;
mutable bool m_analysisIsOk, m_factorizationIsOk;
};
#ifndef EIGEN_TEST_EVALUATORS
namespace internal {
template<typename Derived, typename Rhs>
@@ -248,6 +276,7 @@ struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs>
};
} // end namespace internal
#endif // EIGEN_TEST_EVALUATORS
} // end namespace Eigen

114
Eigen/src/IterativeLinearSolvers/SolveWithGuess.h Normal file
View File

@@ -0,0 +1,114 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_SOLVEWITHGUESS_H
#define EIGEN_SOLVEWITHGUESS_H
namespace Eigen {
template<typename Decomposition, typename RhsType, typename GuessType> class SolveWithGuess;
/** \class SolveWithGuess
* \ingroup IterativeLinearSolvers_Module
*
* \brief Pseudo expression representing a solving operation
*
* \tparam Decomposition the type of the matrix or decomposion object
* \tparam Rhstype the type of the right-hand side
*
* This class represents an expression of A.solve(B)
* and most of the time this is the only way it is used.
*
*/
namespace internal {
template<typename Decomposition, typename RhsType, typename GuessType>
struct traits<SolveWithGuess<Decomposition, RhsType, GuessType> >
: traits<Solve<Decomposition,RhsType> >
{};
}
template<typename Decomposition, typename RhsType, typename GuessType>
class SolveWithGuess : public internal::generic_xpr_base<SolveWithGuess<Decomposition,RhsType,GuessType>, MatrixXpr, typename internal::traits<RhsType>::StorageKind>::type
{
public:
typedef typename RhsType::Index Index;
typedef typename internal::traits<SolveWithGuess>::PlainObject PlainObject;
typedef typename internal::generic_xpr_base<SolveWithGuess<Decomposition,RhsType,GuessType>, MatrixXpr, typename internal::traits<RhsType>::StorageKind>::type Base;
SolveWithGuess(const Decomposition &dec, const RhsType &rhs, const GuessType &guess)
: m_dec(dec), m_rhs(rhs), m_guess(guess)
{}
EIGEN_DEVICE_FUNC Index rows() const { return m_dec.cols(); }
EIGEN_DEVICE_FUNC Index cols() const { return m_rhs.cols(); }
EIGEN_DEVICE_FUNC const Decomposition& dec() const { return m_dec; }
EIGEN_DEVICE_FUNC const RhsType& rhs() const { return m_rhs; }
EIGEN_DEVICE_FUNC const GuessType& guess() const { return m_guess; }
protected:
const Decomposition &m_dec;
const RhsType &m_rhs;
const GuessType &m_guess;
typedef typename internal::traits<SolveWithGuess>::Scalar Scalar;
private:
Scalar coeff(Index row, Index col) const;
Scalar coeff(Index i) const;
};
namespace internal {
// Evaluator of SolveWithGuess -> eval into a temporary
template<typename Decomposition, typename RhsType, typename GuessType>
struct evaluator<SolveWithGuess<Decomposition,RhsType, GuessType> >
: public evaluator<typename SolveWithGuess<Decomposition,RhsType,GuessType>::PlainObject>::type
{
typedef SolveWithGuess<Decomposition,RhsType,GuessType> SolveType;
typedef typename SolveType::PlainObject PlainObject;
typedef typename evaluator<PlainObject>::type Base;
typedef evaluator type;
typedef evaluator nestedType;
evaluator(const SolveType& solve)
: m_result(solve.rows(), solve.cols())
{
::new (static_cast<Base*>(this)) Base(m_result);
solve.dec()._solve_with_guess_impl(solve.rhs(), m_result, solve().guess());
}
protected:
PlainObject m_result;
};
// Specialization for "dst = dec.solve(rhs)"
// NOTE we need to specialize it for Dense2Dense to avoid ambiguous specialization error and a Sparse2Sparse specialization must exist somewhere
template<typename DstXprType, typename DecType, typename RhsType, typename GuessType, typename Scalar>
struct Assignment<DstXprType, SolveWithGuess<DecType,RhsType,GuessType>, internal::assign_op<Scalar>, Dense2Dense, Scalar>
{
typedef Solve<DecType,RhsType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
{
// FIXME shall we resize dst here?
dst = src.guess();
src.dec()._solve_with_guess_impl(src.rhs(), dst/*, src.guess()*/);
}
};
} // end namepsace internal
} // end namespace Eigen
#endif // EIGEN_SOLVEWITHGUESS_H