Add eigen to thirdparty
This commit is contained in:
kradchen
697
thirdparty/include/Eigen/src/SparseCholesky/SimplicialCholesky.h
vendored
Normal file
697
thirdparty/include/Eigen/src/SparseCholesky/SimplicialCholesky.h
vendored
Normal file
@@ -0,0 +1,697 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2008-2012 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_SIMPLICIAL_CHOLESKY_H
|
||||
#define EIGEN_SIMPLICIAL_CHOLESKY_H
|
||||
|
||||
namespace Eigen {
|
||||
|
||||
enum SimplicialCholeskyMode {
|
||||
SimplicialCholeskyLLT,
|
||||
SimplicialCholeskyLDLT
|
||||
};
|
||||
|
||||
namespace internal {
|
||||
template<typename CholMatrixType, typename InputMatrixType>
|
||||
struct simplicial_cholesky_grab_input {
|
||||
typedef CholMatrixType const * ConstCholMatrixPtr;
|
||||
static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp)
|
||||
{
|
||||
tmp = input;
|
||||
pmat = &tmp;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename MatrixType>
|
||||
struct simplicial_cholesky_grab_input<MatrixType,MatrixType> {
|
||||
typedef MatrixType const * ConstMatrixPtr;
|
||||
static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/)
|
||||
{
|
||||
pmat = &input;
|
||||
}
|
||||
};
|
||||
} // end namespace internal
|
||||
|
||||
/** \ingroup SparseCholesky_Module
|
||||
* \brief A base class for direct sparse Cholesky factorizations
|
||||
*
|
||||
* This is a base class for LL^T and LDL^T Cholesky factorizations of sparse matrices that are
|
||||
* selfadjoint and positive definite. These factorizations allow for solving A.X = B where
|
||||
* X and B can be either dense or sparse.
|
||||
*
|
||||
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
|
||||
* such that the factorized matrix is P A P^-1.
|
||||
*
|
||||
* \tparam Derived the type of the derived class, that is the actual factorization type.
|
||||
*
|
||||
*/
|
||||
template<typename Derived>
|
||||
class SimplicialCholeskyBase : public SparseSolverBase<Derived>
|
||||
{
|
||||
typedef SparseSolverBase<Derived> Base;
|
||||
using Base::m_isInitialized;
|
||||
|
||||
public:
|
||||
typedef typename internal::traits<Derived>::MatrixType MatrixType;
|
||||
typedef typename internal::traits<Derived>::OrderingType OrderingType;
|
||||
enum { UpLo = internal::traits<Derived>::UpLo };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
|
||||
typedef CholMatrixType const * ConstCholMatrixPtr;
|
||||
typedef Matrix<Scalar,Dynamic,1> VectorType;
|
||||
typedef Matrix<StorageIndex,Dynamic,1> VectorI;
|
||||
|
||||
enum {
|
||||
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
|
||||
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
|
||||
};
|
||||
|
||||
public:
|
||||
|
||||
using Base::derived;
|
||||
|
||||
/** Default constructor */
|
||||
SimplicialCholeskyBase()
|
||||
: m_info(Success),
|
||||
m_factorizationIsOk(false),
|
||||
m_analysisIsOk(false),
|
||||
m_shiftOffset(0),
|
||||
m_shiftScale(1)
|
||||
{}
|
||||
|
||||
explicit SimplicialCholeskyBase(const MatrixType& matrix)
|
||||
: m_info(Success),
|
||||
m_factorizationIsOk(false),
|
||||
m_analysisIsOk(false),
|
||||
m_shiftOffset(0),
|
||||
m_shiftScale(1)
|
||||
{
|
||||
derived().compute(matrix);
|
||||
}
|
||||
|
||||
~SimplicialCholeskyBase()
|
||||
{
|
||||
}
|
||||
|
||||
Derived& derived() { return *static_cast<Derived*>(this); }
|
||||
const Derived& derived() const { return *static_cast<const Derived*>(this); }
|
||||
|
||||
inline Index cols() const { return m_matrix.cols(); }
|
||||
inline Index rows() const { return m_matrix.rows(); }
|
||||
|
||||
/** \brief Reports whether previous computation was successful.
|
||||
*
|
||||
* \returns \c Success if computation was successful,
|
||||
* \c NumericalIssue if the matrix.appears to be negative.
|
||||
*/
|
||||
ComputationInfo info() const
|
||||
{
|
||||
eigen_assert(m_isInitialized && "Decomposition is not initialized.");
|
||||
return m_info;
|
||||
}
|
||||
|
||||
/** \returns the permutation P
|
||||
* \sa permutationPinv() */
|
||||
const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationP() const
|
||||
{ return m_P; }
|
||||
|
||||
/** \returns the inverse P^-1 of the permutation P
|
||||
* \sa permutationP() */
|
||||
const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationPinv() const
|
||||
{ return m_Pinv; }
|
||||
|
||||
/** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization.
|
||||
*
|
||||
* During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n
|
||||
* \c d_ii = \a offset + \a scale * \c d_ii
|
||||
*
|
||||
* The default is the identity transformation with \a offset=0, and \a scale=1.
|
||||
*
|
||||
* \returns a reference to \c *this.
|
||||
*/
|
||||
Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
|
||||
{
|
||||
m_shiftOffset = offset;
|
||||
m_shiftScale = scale;
|
||||
return derived();
|
||||
}
|
||||
|
||||
#ifndef EIGEN_PARSED_BY_DOXYGEN
|
||||
/** \internal */
|
||||
template<typename Stream>
|
||||
void dumpMemory(Stream& s)
|
||||
{
|
||||
int total = 0;
|
||||
s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
|
||||
s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
|
||||
s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
|
||||
s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
|
||||
s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
|
||||
s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
|
||||
s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
|
||||
}
|
||||
|
||||
/** \internal */
|
||||
template<typename Rhs,typename Dest>
|
||||
void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
|
||||
{
|
||||
eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
|
||||
eigen_assert(m_matrix.rows()==b.rows());
|
||||
|
||||
if(m_info!=Success)
|
||||
return;
|
||||
|
||||
if(m_P.size()>0)
|
||||
dest = m_P * b;
|
||||
else
|
||||
dest = b;
|
||||
|
||||
if(m_matrix.nonZeros()>0) // otherwise L==I
|
||||
derived().matrixL().solveInPlace(dest);
|
||||
|
||||
if(m_diag.size()>0)
|
||||
dest = m_diag.asDiagonal().inverse() * dest;
|
||||
|
||||
if (m_matrix.nonZeros()>0) // otherwise U==I
|
||||
derived().matrixU().solveInPlace(dest);
|
||||
|
||||
if(m_P.size()>0)
|
||||
dest = m_Pinv * dest;
|
||||
}
|
||||
|
||||
template<typename Rhs,typename Dest>
|
||||
void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
|
||||
{
|
||||
internal::solve_sparse_through_dense_panels(derived(), b, dest);
|
||||
}
|
||||
|
||||
#endif // EIGEN_PARSED_BY_DOXYGEN
|
||||
|
||||
protected:
|
||||
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
template<bool DoLDLT>
|
||||
void compute(const MatrixType& matrix)
|
||||
{
|
||||
eigen_assert(matrix.rows()==matrix.cols());
|
||||
Index size = matrix.cols();
|
||||
CholMatrixType tmp(size,size);
|
||||
ConstCholMatrixPtr pmat;
|
||||
ordering(matrix, pmat, tmp);
|
||||
analyzePattern_preordered(*pmat, DoLDLT);
|
||||
factorize_preordered<DoLDLT>(*pmat);
|
||||
}
|
||||
|
||||
template<bool DoLDLT>
|
||||
void factorize(const MatrixType& a)
|
||||
{
|
||||
eigen_assert(a.rows()==a.cols());
|
||||
Index size = a.cols();
|
||||
CholMatrixType tmp(size,size);
|
||||
ConstCholMatrixPtr pmat;
|
||||
|
||||
if(m_P.size() == 0 && (int(UpLo) & int(Upper)) == Upper)
|
||||
{
|
||||
// If there is no ordering, try to directly use the input matrix without any copy
|
||||
internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
|
||||
}
|
||||
else
|
||||
{
|
||||
tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
|
||||
pmat = &tmp;
|
||||
}
|
||||
|
||||
factorize_preordered<DoLDLT>(*pmat);
|
||||
}
|
||||
|
||||
template<bool DoLDLT>
|
||||
void factorize_preordered(const CholMatrixType& a);
|
||||
|
||||
void analyzePattern(const MatrixType& a, bool doLDLT)
|
||||
{
|
||||
eigen_assert(a.rows()==a.cols());
|
||||
Index size = a.cols();
|
||||
CholMatrixType tmp(size,size);
|
||||
ConstCholMatrixPtr pmat;
|
||||
ordering(a, pmat, tmp);
|
||||
analyzePattern_preordered(*pmat,doLDLT);
|
||||
}
|
||||
void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
|
||||
|
||||
void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
|
||||
|
||||
/** keeps off-diagonal entries; drops diagonal entries */
|
||||
struct keep_diag {
|
||||
inline bool operator() (const Index& row, const Index& col, const Scalar&) const
|
||||
{
|
||||
return row!=col;
|
||||
}
|
||||
};
|
||||
|
||||
mutable ComputationInfo m_info;
|
||||
bool m_factorizationIsOk;
|
||||
bool m_analysisIsOk;
|
||||
|
||||
CholMatrixType m_matrix;
|
||||
VectorType m_diag; // the diagonal coefficients (LDLT mode)
|
||||
VectorI m_parent; // elimination tree
|
||||
VectorI m_nonZerosPerCol;
|
||||
PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P; // the permutation
|
||||
PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation
|
||||
|
||||
RealScalar m_shiftOffset;
|
||||
RealScalar m_shiftScale;
|
||||
};
|
||||
|
||||
template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT;
|
||||
template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT;
|
||||
template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky;
|
||||
|
||||
namespace internal {
|
||||
|
||||
template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
typedef _MatrixType MatrixType;
|
||||
typedef _Ordering OrderingType;
|
||||
enum { UpLo = _UpLo };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL;
|
||||
typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
|
||||
static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
|
||||
static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
|
||||
};
|
||||
|
||||
template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
typedef _MatrixType MatrixType;
|
||||
typedef _Ordering OrderingType;
|
||||
enum { UpLo = _UpLo };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
|
||||
typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
|
||||
typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
|
||||
static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
|
||||
static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
|
||||
};
|
||||
|
||||
template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
typedef _MatrixType MatrixType;
|
||||
typedef _Ordering OrderingType;
|
||||
enum { UpLo = _UpLo };
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
/** \ingroup SparseCholesky_Module
|
||||
* \class SimplicialLLT
|
||||
* \brief A direct sparse LLT Cholesky factorizations
|
||||
*
|
||||
* This class provides a LL^T Cholesky factorizations of sparse matrices that are
|
||||
* selfadjoint and positive definite. The factorization allows for solving A.X = B where
|
||||
* X and B can be either dense or sparse.
|
||||
*
|
||||
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
|
||||
* such that the factorized matrix is P A P^-1.
|
||||
*
|
||||
* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
||||
* \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
||||
* or Upper. Default is Lower.
|
||||
* \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
|
||||
*
|
||||
* \implsparsesolverconcept
|
||||
*
|
||||
* \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering
|
||||
*/
|
||||
template<typename _MatrixType, int _UpLo, typename _Ordering>
|
||||
class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
public:
|
||||
typedef _MatrixType MatrixType;
|
||||
enum { UpLo = _UpLo };
|
||||
typedef SimplicialCholeskyBase<SimplicialLLT> Base;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
|
||||
typedef Matrix<Scalar,Dynamic,1> VectorType;
|
||||
typedef internal::traits<SimplicialLLT> Traits;
|
||||
typedef typename Traits::MatrixL MatrixL;
|
||||
typedef typename Traits::MatrixU MatrixU;
|
||||
public:
|
||||
/** Default constructor */
|
||||
SimplicialLLT() : Base() {}
|
||||
/** Constructs and performs the LLT factorization of \a matrix */
|
||||
explicit SimplicialLLT(const MatrixType& matrix)
|
||||
: Base(matrix) {}
|
||||
|
||||
/** \returns an expression of the factor L */
|
||||
inline const MatrixL matrixL() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
|
||||
return Traits::getL(Base::m_matrix);
|
||||
}
|
||||
|
||||
/** \returns an expression of the factor U (= L^*) */
|
||||
inline const MatrixU matrixU() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
|
||||
return Traits::getU(Base::m_matrix);
|
||||
}
|
||||
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
SimplicialLLT& compute(const MatrixType& matrix)
|
||||
{
|
||||
Base::template compute<false>(matrix);
|
||||
return *this;
|
||||
}
|
||||
|
||||
/** Performs a symbolic decomposition on the sparcity of \a matrix.
|
||||
*
|
||||
* This function is particularly useful when solving for several problems having the same structure.
|
||||
*
|
||||
* \sa factorize()
|
||||
*/
|
||||
void analyzePattern(const MatrixType& a)
|
||||
{
|
||||
Base::analyzePattern(a, false);
|
||||
}
|
||||
|
||||
/** Performs a numeric decomposition of \a matrix
|
||||
*
|
||||
* The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
|
||||
*
|
||||
* \sa analyzePattern()
|
||||
*/
|
||||
void factorize(const MatrixType& a)
|
||||
{
|
||||
Base::template factorize<false>(a);
|
||||
}
|
||||
|
||||
/** \returns the determinant of the underlying matrix from the current factorization */
|
||||
Scalar determinant() const
|
||||
{
|
||||
Scalar detL = Base::m_matrix.diagonal().prod();
|
||||
return numext::abs2(detL);
|
||||
}
|
||||
};
|
||||
|
||||
/** \ingroup SparseCholesky_Module
|
||||
* \class SimplicialLDLT
|
||||
* \brief A direct sparse LDLT Cholesky factorizations without square root.
|
||||
*
|
||||
* This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are
|
||||
* selfadjoint and positive definite. The factorization allows for solving A.X = B where
|
||||
* X and B can be either dense or sparse.
|
||||
*
|
||||
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
|
||||
* such that the factorized matrix is P A P^-1.
|
||||
*
|
||||
* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
||||
* \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
||||
* or Upper. Default is Lower.
|
||||
* \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
|
||||
*
|
||||
* \implsparsesolverconcept
|
||||
*
|
||||
* \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering
|
||||
*/
|
||||
template<typename _MatrixType, int _UpLo, typename _Ordering>
|
||||
class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
public:
|
||||
typedef _MatrixType MatrixType;
|
||||
enum { UpLo = _UpLo };
|
||||
typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
|
||||
typedef Matrix<Scalar,Dynamic,1> VectorType;
|
||||
typedef internal::traits<SimplicialLDLT> Traits;
|
||||
typedef typename Traits::MatrixL MatrixL;
|
||||
typedef typename Traits::MatrixU MatrixU;
|
||||
public:
|
||||
/** Default constructor */
|
||||
SimplicialLDLT() : Base() {}
|
||||
|
||||
/** Constructs and performs the LLT factorization of \a matrix */
|
||||
explicit SimplicialLDLT(const MatrixType& matrix)
|
||||
: Base(matrix) {}
|
||||
|
||||
/** \returns a vector expression of the diagonal D */
|
||||
inline const VectorType vectorD() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
||||
return Base::m_diag;
|
||||
}
|
||||
/** \returns an expression of the factor L */
|
||||
inline const MatrixL matrixL() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
||||
return Traits::getL(Base::m_matrix);
|
||||
}
|
||||
|
||||
/** \returns an expression of the factor U (= L^*) */
|
||||
inline const MatrixU matrixU() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
||||
return Traits::getU(Base::m_matrix);
|
||||
}
|
||||
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
SimplicialLDLT& compute(const MatrixType& matrix)
|
||||
{
|
||||
Base::template compute<true>(matrix);
|
||||
return *this;
|
||||
}
|
||||
|
||||
/** Performs a symbolic decomposition on the sparcity of \a matrix.
|
||||
*
|
||||
* This function is particularly useful when solving for several problems having the same structure.
|
||||
*
|
||||
* \sa factorize()
|
||||
*/
|
||||
void analyzePattern(const MatrixType& a)
|
||||
{
|
||||
Base::analyzePattern(a, true);
|
||||
}
|
||||
|
||||
/** Performs a numeric decomposition of \a matrix
|
||||
*
|
||||
* The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
|
||||
*
|
||||
* \sa analyzePattern()
|
||||
*/
|
||||
void factorize(const MatrixType& a)
|
||||
{
|
||||
Base::template factorize<true>(a);
|
||||
}
|
||||
|
||||
/** \returns the determinant of the underlying matrix from the current factorization */
|
||||
Scalar determinant() const
|
||||
{
|
||||
return Base::m_diag.prod();
|
||||
}
|
||||
};
|
||||
|
||||
/** \deprecated use SimplicialLDLT or class SimplicialLLT
|
||||
* \ingroup SparseCholesky_Module
|
||||
* \class SimplicialCholesky
|
||||
*
|
||||
* \sa class SimplicialLDLT, class SimplicialLLT
|
||||
*/
|
||||
template<typename _MatrixType, int _UpLo, typename _Ordering>
|
||||
class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
public:
|
||||
typedef _MatrixType MatrixType;
|
||||
enum { UpLo = _UpLo };
|
||||
typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef typename MatrixType::StorageIndex StorageIndex;
|
||||
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
|
||||
typedef Matrix<Scalar,Dynamic,1> VectorType;
|
||||
typedef internal::traits<SimplicialCholesky> Traits;
|
||||
typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
|
||||
typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
|
||||
public:
|
||||
SimplicialCholesky() : Base(), m_LDLT(true) {}
|
||||
|
||||
explicit SimplicialCholesky(const MatrixType& matrix)
|
||||
: Base(), m_LDLT(true)
|
||||
{
|
||||
compute(matrix);
|
||||
}
|
||||
|
||||
SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
|
||||
{
|
||||
switch(mode)
|
||||
{
|
||||
case SimplicialCholeskyLLT:
|
||||
m_LDLT = false;
|
||||
break;
|
||||
case SimplicialCholeskyLDLT:
|
||||
m_LDLT = true;
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
inline const VectorType vectorD() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
|
||||
return Base::m_diag;
|
||||
}
|
||||
inline const CholMatrixType rawMatrix() const {
|
||||
eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
|
||||
return Base::m_matrix;
|
||||
}
|
||||
|
||||
/** Computes the sparse Cholesky decomposition of \a matrix */
|
||||
SimplicialCholesky& compute(const MatrixType& matrix)
|
||||
{
|
||||
if(m_LDLT)
|
||||
Base::template compute<true>(matrix);
|
||||
else
|
||||
Base::template compute<false>(matrix);
|
||||
return *this;
|
||||
}
|
||||
|
||||
/** Performs a symbolic decomposition on the sparcity of \a matrix.
|
||||
*
|
||||
* This function is particularly useful when solving for several problems having the same structure.
|
||||
*
|
||||
* \sa factorize()
|
||||
*/
|
||||
void analyzePattern(const MatrixType& a)
|
||||
{
|
||||
Base::analyzePattern(a, m_LDLT);
|
||||
}
|
||||
|
||||
/** Performs a numeric decomposition of \a matrix
|
||||
*
|
||||
* The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
|
||||
*
|
||||
* \sa analyzePattern()
|
||||
*/
|
||||
void factorize(const MatrixType& a)
|
||||
{
|
||||
if(m_LDLT)
|
||||
Base::template factorize<true>(a);
|
||||
else
|
||||
Base::template factorize<false>(a);
|
||||
}
|
||||
|
||||
/** \internal */
|
||||
template<typename Rhs,typename Dest>
|
||||
void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
|
||||
{
|
||||
eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
|
||||
eigen_assert(Base::m_matrix.rows()==b.rows());
|
||||
|
||||
if(Base::m_info!=Success)
|
||||
return;
|
||||
|
||||
if(Base::m_P.size()>0)
|
||||
dest = Base::m_P * b;
|
||||
else
|
||||
dest = b;
|
||||
|
||||
if(Base::m_matrix.nonZeros()>0) // otherwise L==I
|
||||
{
|
||||
if(m_LDLT)
|
||||
LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
|
||||
else
|
||||
LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
|
||||
}
|
||||
|
||||
if(Base::m_diag.size()>0)
|
||||
dest = Base::m_diag.real().asDiagonal().inverse() * dest;
|
||||
|
||||
if (Base::m_matrix.nonZeros()>0) // otherwise I==I
|
||||
{
|
||||
if(m_LDLT)
|
||||
LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
|
||||
else
|
||||
LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
|
||||
}
|
||||
|
||||
if(Base::m_P.size()>0)
|
||||
dest = Base::m_Pinv * dest;
|
||||
}
|
||||
|
||||
/** \internal */
|
||||
template<typename Rhs,typename Dest>
|
||||
void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
|
||||
{
|
||||
internal::solve_sparse_through_dense_panels(*this, b, dest);
|
||||
}
|
||||
|
||||
Scalar determinant() const
|
||||
{
|
||||
if(m_LDLT)
|
||||
{
|
||||
return Base::m_diag.prod();
|
||||
}
|
||||
else
|
||||
{
|
||||
Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
|
||||
return numext::abs2(detL);
|
||||
}
|
||||
}
|
||||
|
||||
protected:
|
||||
bool m_LDLT;
|
||||
};
|
||||
|
||||
template<typename Derived>
|
||||
void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
|
||||
{
|
||||
eigen_assert(a.rows()==a.cols());
|
||||
const Index size = a.rows();
|
||||
pmat = ≈
|
||||
// Note that ordering methods compute the inverse permutation
|
||||
if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
|
||||
{
|
||||
{
|
||||
CholMatrixType C;
|
||||
C = a.template selfadjointView<UpLo>();
|
||||
|
||||
OrderingType ordering;
|
||||
ordering(C,m_Pinv);
|
||||
}
|
||||
|
||||
if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
|
||||
else m_P.resize(0);
|
||||
|
||||
ap.resize(size,size);
|
||||
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
|
||||
}
|
||||
else
|
||||
{
|
||||
m_Pinv.resize(0);
|
||||
m_P.resize(0);
|
||||
if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor)
|
||||
{
|
||||
// we have to transpose the lower part to to the upper one
|
||||
ap.resize(size,size);
|
||||
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
|
||||
}
|
||||
else
|
||||
internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
|
||||
}
|
||||
}
|
||||
|
||||
} // end namespace Eigen
|
||||
|
||||
#endif // EIGEN_SIMPLICIAL_CHOLESKY_H
|
||||
174
thirdparty/include/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h
vendored
Normal file
174
thirdparty/include/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h
vendored
Normal file
@@ -0,0 +1,174 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2008-2012 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/.
|
||||
|
||||
/*
|
||||
NOTE: these functions have been adapted from the LDL library:
|
||||
|
||||
LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved.
|
||||
|
||||
The author of LDL, Timothy A. Davis., has executed a license with Google LLC
|
||||
to permit distribution of this code and derivative works as part of Eigen under
|
||||
the Mozilla Public License v. 2.0, as stated at the top of this file.
|
||||
*/
|
||||
|
||||
#ifndef EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H
|
||||
#define EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H
|
||||
|
||||
namespace Eigen {
|
||||
|
||||
template<typename Derived>
|
||||
void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT)
|
||||
{
|
||||
const StorageIndex size = StorageIndex(ap.rows());
|
||||
m_matrix.resize(size, size);
|
||||
m_parent.resize(size);
|
||||
m_nonZerosPerCol.resize(size);
|
||||
|
||||
ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0);
|
||||
|
||||
for(StorageIndex k = 0; k < size; ++k)
|
||||
{
|
||||
/* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
|
||||
m_parent[k] = -1; /* parent of k is not yet known */
|
||||
tags[k] = k; /* mark node k as visited */
|
||||
m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */
|
||||
for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
|
||||
{
|
||||
StorageIndex i = it.index();
|
||||
if(i < k)
|
||||
{
|
||||
/* follow path from i to root of etree, stop at flagged node */
|
||||
for(; tags[i] != k; i = m_parent[i])
|
||||
{
|
||||
/* find parent of i if not yet determined */
|
||||
if (m_parent[i] == -1)
|
||||
m_parent[i] = k;
|
||||
m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */
|
||||
tags[i] = k; /* mark i as visited */
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* construct Lp index array from m_nonZerosPerCol column counts */
|
||||
StorageIndex* Lp = m_matrix.outerIndexPtr();
|
||||
Lp[0] = 0;
|
||||
for(StorageIndex k = 0; k < size; ++k)
|
||||
Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
|
||||
|
||||
m_matrix.resizeNonZeros(Lp[size]);
|
||||
|
||||
m_isInitialized = true;
|
||||
m_info = Success;
|
||||
m_analysisIsOk = true;
|
||||
m_factorizationIsOk = false;
|
||||
}
|
||||
|
||||
|
||||
template<typename Derived>
|
||||
template<bool DoLDLT>
|
||||
void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap)
|
||||
{
|
||||
using std::sqrt;
|
||||
|
||||
eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
|
||||
eigen_assert(ap.rows()==ap.cols());
|
||||
eigen_assert(m_parent.size()==ap.rows());
|
||||
eigen_assert(m_nonZerosPerCol.size()==ap.rows());
|
||||
|
||||
const StorageIndex size = StorageIndex(ap.rows());
|
||||
const StorageIndex* Lp = m_matrix.outerIndexPtr();
|
||||
StorageIndex* Li = m_matrix.innerIndexPtr();
|
||||
Scalar* Lx = m_matrix.valuePtr();
|
||||
|
||||
ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
|
||||
ei_declare_aligned_stack_constructed_variable(StorageIndex, pattern, size, 0);
|
||||
ei_declare_aligned_stack_constructed_variable(StorageIndex, tags, size, 0);
|
||||
|
||||
bool ok = true;
|
||||
m_diag.resize(DoLDLT ? size : 0);
|
||||
|
||||
for(StorageIndex k = 0; k < size; ++k)
|
||||
{
|
||||
// compute nonzero pattern of kth row of L, in topological order
|
||||
y[k] = Scalar(0); // Y(0:k) is now all zero
|
||||
StorageIndex top = size; // stack for pattern is empty
|
||||
tags[k] = k; // mark node k as visited
|
||||
m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L
|
||||
for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
|
||||
{
|
||||
StorageIndex i = it.index();
|
||||
if(i <= k)
|
||||
{
|
||||
y[i] += numext::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */
|
||||
Index len;
|
||||
for(len = 0; tags[i] != k; i = m_parent[i])
|
||||
{
|
||||
pattern[len++] = i; /* L(k,i) is nonzero */
|
||||
tags[i] = k; /* mark i as visited */
|
||||
}
|
||||
while(len > 0)
|
||||
pattern[--top] = pattern[--len];
|
||||
}
|
||||
}
|
||||
|
||||
/* compute numerical values kth row of L (a sparse triangular solve) */
|
||||
|
||||
RealScalar d = numext::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k)
|
||||
y[k] = Scalar(0);
|
||||
for(; top < size; ++top)
|
||||
{
|
||||
Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */
|
||||
Scalar yi = y[i]; /* get and clear Y(i) */
|
||||
y[i] = Scalar(0);
|
||||
|
||||
/* the nonzero entry L(k,i) */
|
||||
Scalar l_ki;
|
||||
if(DoLDLT)
|
||||
l_ki = yi / numext::real(m_diag[i]);
|
||||
else
|
||||
yi = l_ki = yi / Lx[Lp[i]];
|
||||
|
||||
Index p2 = Lp[i] + m_nonZerosPerCol[i];
|
||||
Index p;
|
||||
for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
|
||||
y[Li[p]] -= numext::conj(Lx[p]) * yi;
|
||||
d -= numext::real(l_ki * numext::conj(yi));
|
||||
Li[p] = k; /* store L(k,i) in column form of L */
|
||||
Lx[p] = l_ki;
|
||||
++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
|
||||
}
|
||||
if(DoLDLT)
|
||||
{
|
||||
m_diag[k] = d;
|
||||
if(d == RealScalar(0))
|
||||
{
|
||||
ok = false; /* failure, D(k,k) is zero */
|
||||
break;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
Index p = Lp[k] + m_nonZerosPerCol[k]++;
|
||||
Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */
|
||||
if(d <= RealScalar(0)) {
|
||||
ok = false; /* failure, matrix is not positive definite */
|
||||
break;
|
||||
}
|
||||
Lx[p] = sqrt(d) ;
|
||||
}
|
||||
}
|
||||
|
||||
m_info = ok ? Success : NumericalIssue;
|
||||
m_factorizationIsOk = true;
|
||||
}
|
||||
|
||||
} // end namespace Eigen
|
||||
|
||||
#endif // EIGEN_SIMPLICIAL_CHOLESKY_IMPL_H
|
||||
Reference in New Issue
Block a user