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

fix sparse rankUpdate and triangularView iterator

Browse Source
This commit is contained in:
Gael Guennebaud 2010-10-27 15:13:03 +02:00
parent 241e5ee3e7
commit 02c8b6af82
3 changed files with 24 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Eigen/src/Sparse/SparseSelfAdjointView.h
View File

@ -92,7 +92,7 @@ template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
* call this function with u.adjoint(). * call this function with u.adjoint().
*/ */
template<typename DerivedU> template<typename DerivedU>
SparseSelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1)); SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
// const SparseLLT<PlainObject, UpLo> llt() const; // const SparseLLT<PlainObject, UpLo> llt() const;
// const SparseLDLT<PlainObject, UpLo> ldlt() const; // const SparseLDLT<PlainObject, UpLo> ldlt() const;
@ -127,15 +127,15 @@ SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView(
template<typename MatrixType, unsigned int UpLo> template<typename MatrixType, unsigned int UpLo>
template<typename DerivedU> template<typename DerivedU>
SparseSelfAdjointView<MatrixType,UpLo>& SparseSelfAdjointView<MatrixType,UpLo>&
SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha) SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha)
{ {
SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint(); SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint();
if(alpha==Scalar(0)) if(alpha==Scalar(0))
m_matrix = tmp.template triangularView<UpLo>(); m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
else else
m_matrix += alpha * tmp.template triangularView<UpLo>(); m_matrix.const_cast_derived() /*+*/= alpha * tmp.template triangularView<UpLo>();
return this; return *this;
} }
/*************************************************************************** /***************************************************************************

16
Eigen/src/Sparse/SparseTriangularView.h
View File

@ -26,11 +26,13 @@
#define EIGEN_SPARSE_TRIANGULARVIEW_H #define EIGEN_SPARSE_TRIANGULARVIEW_H
namespace internal { namespace internal {
template<typename MatrixType, int Mode> template<typename MatrixType, int Mode>
struct traits<SparseTriangularView<MatrixType,Mode> > struct traits<SparseTriangularView<MatrixType,Mode> >
: public traits<MatrixType> : public traits<MatrixType>
{}; {};
}
} // namespace internal
template<typename MatrixType, int Mode> class SparseTriangularView template<typename MatrixType, int Mode> class SparseTriangularView
: public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> > : public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> >
@ -39,12 +41,12 @@ template<typename MatrixType, int Mode> class SparseTriangularView
|| (Mode==Upper && (MatrixType::Flags&RowMajorBit)) }; || (Mode==Upper && (MatrixType::Flags&RowMajorBit)) };
public: public:
class InnerIterator; EIGEN_SPARSE_PUBLIC_INTERFACE(SparseTriangularView)
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index;
inline Index rows() { return m_matrix.rows(); } class InnerIterator;
inline Index cols() { return m_matrix.cols(); }
inline Index rows() const { return m_matrix.rows(); }
inline Index cols() const { return m_matrix.cols(); }
typedef typename internal::conditional<internal::must_nest_by_value<MatrixType>::ret, typedef typename internal::conditional<internal::must_nest_by_value<MatrixType>::ret,
MatrixType, const MatrixType&>::type MatrixTypeNested; MatrixType, const MatrixType&>::type MatrixTypeNested;
@ -83,7 +85,7 @@ class SparseTriangularView<MatrixType,Mode>::InnerIterator : public MatrixType::
EIGEN_STRONG_INLINE operator bool() const EIGEN_STRONG_INLINE operator bool() const
{ {
return SkipFirst ? Base::operator bool() : (Base::operator bool() && this->index() < this->outer()); return SkipFirst ? Base::operator bool() : (Base::operator bool() && this->index() <= this->outer());
} }
}; };

13
unsupported/test/sparse_llt.cpp
View File

@ -77,8 +77,11 @@ template<typename Scalar> void sparse_llt(int rows, int cols)
// new API // new API
{ {
// Cholmod, as configured in CholmodSupport.h, only supports self-adjoint matrices // Cholmod, as configured in CholmodSupport.h, only supports self-adjoint matrices
SparseMatrix<Scalar> m3 = m2.adjoint()*m2; SparseMatrix<Scalar> m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
DenseMatrix refMat3 = refMat2.adjoint()*refMat2; DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
m3_lo.template selfadjointView<Lower>().rankUpdate(m2,0);
m3_up.template selfadjointView<Upper>().rankUpdate(m2,0);
// with a single vector as the rhs // with a single vector as the rhs
ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b); ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b);
@ -89,6 +92,12 @@ template<typename Scalar> void sparse_llt(int rows, int cols)
x = CholmodDecomposition<SparseMatrix<Scalar>, Upper>(m3).solve(b); x = CholmodDecomposition<SparseMatrix<Scalar>, Upper>(m3).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs"); VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
x = CholmodDecomposition<SparseMatrix<Scalar>, Lower>(m3_lo).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
x = CholmodDecomposition<SparseMatrix<Scalar>, Upper>(m3_up).solve(b);
VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
// with multiple rhs // with multiple rhs
ref_X = refMat3.template selfadjointView<Lower>().llt().solve(B); ref_X = refMat3.template selfadjointView<Lower>().llt().solve(B);