diff options
Diffstat (limited to 'Eigen/src')
| -rw-r--r-- | Eigen/src/Core/MatrixBase.h | 9 | ||||
| -rw-r--r-- | Eigen/src/Core/SelfAdjointView.h | 8 | ||||
| -rw-r--r-- | Eigen/src/Core/util/Constants.h | 2 | ||||
| -rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 4 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseCwiseBinaryOp.h | 36 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseLDLT.h | 4 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseLLT.h | 6 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseMatrix.h | 1 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseMatrixBase.h | 14 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseProduct.h | 186 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseSelfAdjointView.h | 223 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseTriangular.h | 60 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseTriangularView.h | 95 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseUtil.h | 30 | ||||
| -rw-r--r-- | Eigen/src/Sparse/TriangularSolver.h | 18 |
15 files changed, 437 insertions, 259 deletions
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index b752c7821..67e270af7 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -710,15 +710,6 @@ template<typename Derived> class MatrixBase typedef Homogeneous<Derived,MatrixBase<Derived>::ColsAtCompileTime==1?Vertical:Horizontal> HomogeneousReturnType; const HomogeneousReturnType homogeneous() const; -/////////// Sparse module /////////// - - // dense = sparse * dense - template<typename Derived1, typename Derived2> - Derived& lazyAssign(const SparseProduct<Derived1,Derived2,SparseTimeDenseProduct>& product); - // dense = dense * sparse - template<typename Derived1, typename Derived2> - Derived& lazyAssign(const SparseProduct<Derived1,Derived2,DenseTimeSparseProduct>& product); - ////////// Householder module /////////// void makeHouseholderInPlace(Scalar& tau, RealScalar& beta); diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 95ce666c9..3e435853c 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -201,15 +201,15 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, Dynami ***************************************************************************/ template<typename Derived> -template<unsigned int Mode> -const SelfAdjointView<Derived, Mode> MatrixBase<Derived>::selfadjointView() const +template<unsigned int UpLo> +const SelfAdjointView<Derived, UpLo> MatrixBase<Derived>::selfadjointView() const { return derived(); } template<typename Derived> -template<unsigned int Mode> -SelfAdjointView<Derived, Mode> MatrixBase<Derived>::selfadjointView() +template<unsigned int UpLo> +SelfAdjointView<Derived, UpLo> MatrixBase<Derived>::selfadjointView() { return derived(); } diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index ed1cfa35c..c801d8049 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -193,7 +193,7 @@ enum { AsRequested=0, EnforceAlignedAccess=2 }; enum { ConditionalJumpCost = 5 }; enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight }; enum DirectionType { Vertical, Horizontal, BothDirections }; -enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct }; +enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct }; enum { /** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 5591c754d..2fedbbc07 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -146,8 +146,4 @@ template<typename Scalar,int Dim> class Translation; template<typename Scalar> class UniformScaling; template<typename MatrixType,int Direction> class Homogeneous; -// Sparse module: -template<typename Lhs, typename Rhs, int ProductMode> class SparseProduct; - - #endif // EIGEN_FORWARDDECLARATIONS_H diff --git a/Eigen/src/Sparse/SparseCwiseBinaryOp.h b/Eigen/src/Sparse/SparseCwiseBinaryOp.h index 9c354b4d3..5dd7edc00 100644 --- a/Eigen/src/Sparse/SparseCwiseBinaryOp.h +++ b/Eigen/src/Sparse/SparseCwiseBinaryOp.h @@ -42,35 +42,6 @@ // 4 - dense op dense product dense // generic dense -// template<typename BinaryOp, typename Lhs, typename Rhs> -// struct ei_traits<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> > -// { -// typedef typename ei_result_of< -// BinaryOp( -// typename Lhs::Scalar, -// typename Rhs::Scalar -// ) -// >::type Scalar; -// typedef typename ei_promote_storage_type<typename ei_traits<Lhs>::StorageType, -// typename ei_traits<Rhs>::StorageType>::ret StorageType; -// typedef typename Lhs::Nested LhsNested; -// typedef typename Rhs::Nested RhsNested; -// typedef typename ei_unref<LhsNested>::type _LhsNested; -// typedef typename ei_unref<RhsNested>::type _RhsNested; -// enum { -// LhsCoeffReadCost = _LhsNested::CoeffReadCost, -// RhsCoeffReadCost = _RhsNested::CoeffReadCost, -// LhsFlags = _LhsNested::Flags, -// RhsFlags = _RhsNested::Flags, -// RowsAtCompileTime = Lhs::RowsAtCompileTime, -// ColsAtCompileTime = Lhs::ColsAtCompileTime, -// MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime, -// MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime, -// Flags = (int(LhsFlags) | int(RhsFlags)) & HereditaryBits, -// CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost -// }; -// }; - template<> struct ei_promote_storage_type<Dense,Sparse> { typedef Sparse ret; }; @@ -82,16 +53,9 @@ class CwiseBinaryOpImpl<BinaryOp, Lhs, Rhs, Sparse> : public SparseMatrixBase<CwiseBinaryOp<BinaryOp, Lhs, Rhs> > { public: - class InnerIterator; - typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> Derived; EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) -// typedef typename ei_traits<SparseCwiseBinaryOp>::LhsNested LhsNested; -// typedef typename ei_traits<SparseCwiseBinaryOp>::RhsNested RhsNested; -// typedef typename ei_unref<LhsNested>::type _LhsNested; -// typedef typename ei_unref<RhsNested>::type _RhsNested; - }; template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived, diff --git a/Eigen/src/Sparse/SparseLDLT.h b/Eigen/src/Sparse/SparseLDLT.h index 631bff04b..ea5d345f5 100644 --- a/Eigen/src/Sparse/SparseLDLT.h +++ b/Eigen/src/Sparse/SparseLDLT.h @@ -333,12 +333,12 @@ bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const return false; if (m_matrix.nonZeros()>0) // otherwise L==I - m_matrix.template triangular<UnitLowerTriangular>().solveInPlace(b); + m_matrix.template triangularView<UnitLowerTriangular>().solveInPlace(b); b = b.cwise() / m_diag; // FIXME should be .adjoint() but it fails to compile... if (m_matrix.nonZeros()>0) // otherwise L==I - m_matrix.transpose().template triangular<UnitUpperTriangular>().solveInPlace(b); + m_matrix.transpose().template triangularView<UnitUpperTriangular>().solveInPlace(b); return true; } diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h index b2f65b944..6307b2493 100644 --- a/Eigen/src/Sparse/SparseLLT.h +++ b/Eigen/src/Sparse/SparseLLT.h @@ -193,15 +193,15 @@ bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const const int size = m_matrix.rows(); ei_assert(size==b.rows()); - m_matrix.template triangular<LowerTriangular>().solveInPlace(b); + m_matrix.template triangularView<LowerTriangular>().solveInPlace(b); // FIXME should be simply .adjoint() but it fails to compile... if (NumTraits<Scalar>::IsComplex) { CholMatrixType aux = m_matrix.conjugate(); - aux.transpose().template triangular<UpperTriangular>().solveInPlace(b); + aux.transpose().template triangularView<UpperTriangular>().solveInPlace(b); } else - m_matrix.transpose().template triangular<UpperTriangular>().solveInPlace(b); + m_matrix.transpose().template triangularView<UpperTriangular>().solveInPlace(b); return true; } diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index f3915af70..5e89d3f53 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -540,6 +540,7 @@ class SparseMatrix<Scalar,_Options>::InnerIterator inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.m_data.value(m_id)); } inline int index() const { return m_matrix.m_data.index(m_id); } + inline int outer() const { return m_outer; } inline int row() const { return IsRowMajor ? m_outer : index(); } inline int col() const { return IsRowMajor ? index() : m_outer; } diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 2eda425a7..2b6970818 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -244,7 +244,7 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived } template<typename Lhs, typename Rhs> - inline Derived& operator=(const SparseProduct<Lhs,Rhs,SparseTimeSparseProduct>& product); + inline Derived& operator=(const SparseProduct<Lhs,Rhs>& product); friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m) { @@ -337,12 +337,13 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived // dense * sparse (return a dense object) template<typename OtherDerived> friend - const typename SparseProductReturnType<OtherDerived,Derived>::Type + const DenseTimeSparseProduct<OtherDerived,Derived> operator*(const MatrixBase<OtherDerived>& lhs, const Derived& rhs) - { return typename SparseProductReturnType<OtherDerived,Derived>::Type(lhs.derived(),rhs); } + { return DenseTimeSparseProduct<OtherDerived,Derived>(lhs.derived(),rhs); } + // sparse * dense (returns a dense object) template<typename OtherDerived> - const typename SparseProductReturnType<Derived,OtherDerived>::Type + const SparseTimeDenseProduct<Derived,OtherDerived> operator*(const MatrixBase<OtherDerived> &other) const; template<typename OtherDerived> @@ -360,7 +361,10 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived // void solveTriangularInPlace(SparseMatrixBase<OtherDerived>& other) const; template<int Mode> - inline const SparseTriangular<Derived, Mode> triangular() const; + inline const SparseTriangularView<Derived, Mode> triangularView() const; + + template<unsigned int UpLo> inline const SparseSelfAdjointView<Derived, UpLo> selfadjointView() const; + template<unsigned int UpLo> inline SparseSelfAdjointView<Derived, UpLo> selfadjointView(); template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const; template<typename OtherDerived> Scalar dot(const SparseMatrixBase<OtherDerived>& other) const; diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h index 16f36d640..2aaa8713c 100644 --- a/Eigen/src/Sparse/SparseProduct.h +++ b/Eigen/src/Sparse/SparseProduct.h @@ -25,22 +25,8 @@ #ifndef EIGEN_SPARSEPRODUCT_H #define EIGEN_SPARSEPRODUCT_H -template<typename Lhs, typename Rhs> struct ei_sparse_product_mode<Lhs,Rhs,Sparse,Sparse> { enum { value = SparseTimeSparseProduct }; }; -template<typename Lhs, typename Rhs> struct ei_sparse_product_mode<Lhs,Rhs,Dense, Sparse> { enum { value = DenseTimeSparseProduct }; }; -template<typename Lhs, typename Rhs> struct ei_sparse_product_mode<Lhs,Rhs,Sparse,Dense> { enum { value = SparseTimeDenseProduct }; }; - -template<typename Lhs, typename Rhs, int ProductMode> -struct SparseProductReturnType -{ - typedef const typename ei_nested<Lhs,Rhs::RowsAtCompileTime>::type LhsNested; - typedef const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested; - - typedef SparseProduct<LhsNested, RhsNested, ProductMode> Type; -}; - -// sparse product return type specialization template<typename Lhs, typename Rhs> -struct SparseProductReturnType<Lhs,Rhs,SparseTimeSparseProduct> +struct SparseProductReturnType { typedef typename ei_traits<Lhs>::Scalar Scalar; enum { @@ -60,11 +46,11 @@ struct SparseProductReturnType<Lhs,Rhs,SparseTimeSparseProduct> SparseMatrix<Scalar,0>, const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type>::ret RhsNested; - typedef SparseProduct<LhsNested, RhsNested, SparseTimeSparseProduct> Type; + typedef SparseProduct<LhsNested, RhsNested> Type; }; -template<typename LhsNested, typename RhsNested, int ProductMode> -struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> > +template<typename LhsNested, typename RhsNested> +struct ei_traits<SparseProduct<LhsNested, RhsNested> > { // clean the nested types: typedef typename ei_cleantype<LhsNested>::type _LhsNested; @@ -85,7 +71,6 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> > MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime, EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit), - ResultIsSparse = ProductMode==SparseTimeSparseProduct, RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), @@ -96,16 +81,14 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> > CoeffReadCost = Dynamic }; - typedef typename ei_meta_if<ResultIsSparse, Sparse, Dense>::ret StorageType; + typedef Sparse StorageType; - typedef typename ei_meta_if<ResultIsSparse, - SparseMatrixBase<SparseProduct<LhsNested, RhsNested, ProductMode> >, - MatrixBase<SparseProduct<LhsNested, RhsNested, ProductMode> > >::ret Base; + typedef SparseMatrixBase<SparseProduct<LhsNested, RhsNested> > Base; }; -template<typename LhsNested, typename RhsNested, int ProductMode> +template<typename LhsNested, typename RhsNested> class SparseProduct : ei_no_assignment_operator, - public ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >::Base + public ei_traits<SparseProduct<LhsNested, RhsNested> >::Base { public: @@ -256,38 +239,14 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor> } }; -// NOTE eventually let's transpose one argument even in this case since it might be expensive if -// the result is not dense. -// template<typename Lhs, typename Rhs, typename ResultType, int ResStorageOrder> -// struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ResStorageOrder> -// { -// static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) -// { -// // trivial product as lhs.row/rhs.col dot products -// // loop over the preferred order of the result -// } -// }; - // NOTE the 2 others cases (col row *) must never occurs since they are caught // by ProductReturnType which transform it to (col col *) by evaluating rhs. -// template<typename Derived> -// template<typename Lhs, typename Rhs> -// inline Derived& SparseMatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs>& product) -// { -// // std::cout << "sparse product to dense\n"; -// ei_sparse_product_selector< -// typename ei_cleantype<Lhs>::type, -// typename ei_cleantype<Rhs>::type, -// typename ei_cleantype<Derived>::type>::run(product.lhs(),product.rhs(),derived()); -// return derived(); -// } - // sparse = sparse * sparse template<typename Derived> template<typename Lhs, typename Rhs> -inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs,SparseTimeSparseProduct>& product) +inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs>& product) { ei_sparse_product_selector< typename ei_cleantype<Lhs>::type, @@ -296,78 +255,79 @@ inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs return derived(); } -// dense = sparse * dense -// Note that here we force no inlining and separate the setZero() because GCC messes up otherwise -template<typename Lhs, typename Rhs, typename Dest> -EIGEN_DONT_INLINE void ei_sparse_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) + +template<typename Lhs, typename Rhs> +struct ei_traits<SparseTimeDenseProduct<Lhs,Rhs> > + : ei_traits<ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> > { - typedef typename ei_cleantype<Lhs>::type _Lhs; - typedef typename ei_cleantype<Rhs>::type _Rhs; - typedef typename _Lhs::InnerIterator LhsInnerIterator; - enum { - LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, - LhsIsSelfAdjoint = (_Lhs::Flags&SelfAdjointBit)==SelfAdjointBit, - ProcessFirstHalf = LhsIsSelfAdjoint - && ( ((_Lhs::Flags&(UpperTriangularBit|LowerTriangularBit))==0) - || ( (_Lhs::Flags&UpperTriangularBit) && !LhsIsRowMajor) - || ( (_Lhs::Flags&LowerTriangularBit) && LhsIsRowMajor) ), - ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf) - }; - for (int j=0; j<lhs.outerSize(); ++j) - { - LhsInnerIterator i(lhs,j); - if (ProcessSecondHalf && i && (i.index()==j)) - { - dst.row(j) += i.value() * rhs.row(j); - ++i; - } - typename Rhs::Scalar rhs_j = rhs.coeff(j,0); - Block<Dest,1,Dest::ColsAtCompileTime> dst_j(dst.row(LhsIsRowMajor ? j : 0)); - for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) + typedef Dense StorageType; +}; + +template<typename Lhs, typename Rhs> +class SparseTimeDenseProduct + : public ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseTimeDenseProduct) + + SparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const { - if(LhsIsSelfAdjoint) + typedef typename ei_cleantype<Lhs>::type _Lhs; + typedef typename ei_cleantype<Rhs>::type _Rhs; + typedef typename _Lhs::InnerIterator LhsInnerIterator; + enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit }; + for(int j=0; j<m_lhs.outerSize(); ++j) { - int a = LhsIsRowMajor ? j : i.index(); - int b = LhsIsRowMajor ? i.index() : j; - typename Lhs::Scalar v = i.value(); - dst.row(a) += (v) * rhs.row(b); - dst.row(b) += ei_conj(v) * rhs.row(a); + typename Rhs::Scalar rhs_j = alpha * m_rhs.coeff(j,0); + Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0)); + for(LhsInnerIterator it(m_lhs,j); it ;++it) + { + if(LhsIsRowMajor) dest_j += (alpha*it.value()) * m_rhs.row(it.index()); + else if(Rhs::ColsAtCompileTime==1) dest.coeffRef(it.index()) += it.value() * rhs_j; + else dest.row(it.index()) += (alpha*it.value()) * m_rhs.row(j); + } } - else if(LhsIsRowMajor) - dst_j += i.value() * rhs.row(i.index()); - else if(Rhs::ColsAtCompileTime==1) - dst.coeffRef(i.index()) += i.value() * rhs_j; - else - dst.row(i.index()) += i.value() * rhs.row(j); } - if (ProcessFirstHalf && i && (i.index()==j)) - dst.row(j) += i.value() * rhs.row(j); - } -} -template<typename Derived> + private: + SparseTimeDenseProduct& operator=(const SparseTimeDenseProduct&); +}; + + +// dense = dense * sparse template<typename Lhs, typename Rhs> -Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product) +struct ei_traits<DenseTimeSparseProduct<Lhs,Rhs> > + : ei_traits<ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs> > { - derived().setZero(); - ei_sparse_time_dense_product(product.lhs(), product.rhs(), derived()); - return derived(); -} + typedef Dense StorageType; +}; -// dense = dense * sparse -template<typename Derived> template<typename Lhs, typename Rhs> -Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,DenseTimeSparseProduct>& product) +class DenseTimeSparseProduct + : public ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs> { - typedef typename ei_cleantype<Rhs>::type _Rhs; - typedef typename _Rhs::InnerIterator RhsInnerIterator; - enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit }; - derived().setZero(); - for (int j=0; j<product.rhs().outerSize(); ++j) - for (RhsInnerIterator i(product.rhs(),j); i; ++i) - derived().col(RhsIsRowMajor ? i.index() : j) += i.value() * product.lhs().col(RhsIsRowMajor ? j : i.index()); - return derived(); -} + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseProduct) + + DenseTimeSparseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const + { + typedef typename ei_cleantype<Rhs>::type _Rhs; + typedef typename _Rhs::InnerIterator RhsInnerIterator; + enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit }; + for(int j=0; j<m_rhs.outerSize(); ++j) + for(RhsInnerIterator i(m_rhs,j); i; ++i) + dest.col(RhsIsRowMajor ? i.index() : j) += (alpha*i.value()) * m_lhs.col(RhsIsRowMajor ? j : i.index()); + } + + private: + DenseTimeSparseProduct& operator=(const DenseTimeSparseProduct&); +}; // sparse * sparse template<typename Derived> @@ -381,10 +341,10 @@ SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other // sparse * dense template<typename Derived> template<typename OtherDerived> -EIGEN_STRONG_INLINE const typename SparseProductReturnType<Derived,OtherDerived>::Type +EIGEN_STRONG_INLINE const SparseTimeDenseProduct<Derived,OtherDerived> SparseMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const { - return typename SparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); + return SparseTimeDenseProduct<Derived,OtherDerived>(derived(), other.derived()); } #endif // EIGEN_SPARSEPRODUCT_H diff --git a/Eigen/src/Sparse/SparseSelfAdjointView.h b/Eigen/src/Sparse/SparseSelfAdjointView.h new file mode 100644 index 000000000..f5296accf --- /dev/null +++ b/Eigen/src/Sparse/SparseSelfAdjointView.h @@ -0,0 +1,223 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_SELFADJOINTVIEW_H +#define EIGEN_SPARSE_SELFADJOINTVIEW_H + +/** \class SparseSelfAdjointView + * \nonstableyet + * + * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix. + * + * \param MatrixType the type of the dense matrix storing the coefficients + * \param UpLo can be either \c LowerTriangular or \c UpperTriangular + * + * This class is an expression of a sefladjoint matrix from a triangular part of a matrix + * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() + * and most of the time this is the only way that it is used. + * + * \sa SparseMatrixBase::selfAdjointView() + */ +template<typename Lhs, typename Rhs, int UpLo> +class SparseSelfAdjointTimeDenseProduct; + +template<typename Lhs, typename Rhs, int UpLo> +class DenseTimeSparseSelfAdjointProduct; + +template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView +{ + public: + + typedef typename ei_traits<MatrixType>::Scalar Scalar; + + inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix) + { + ei_assert(ei_are_flags_consistent<UpLo>::ret); + ei_assert(rows()==cols() && "SelfAdjointView is only for squared matrices"); + } + + inline int rows() const { return m_matrix.rows(); } + inline int cols() const { return m_matrix.cols(); } + + /** \internal \returns a reference to the nested matrix */ + const MatrixType& matrix() const { return m_matrix; } + + /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ + template<typename OtherDerived> + SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo> + operator*(const MatrixBase<OtherDerived>& rhs) const + { + return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived()); + } + + /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ + template<typename OtherDerived> friend + DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo> + operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs) + { + return DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>(lhs.derived(), rhs.m_matrix); + } + + /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: + * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. + * + * \returns a reference to \c *this + * + * Note that it is faster to set alpha=0 than initializing the matrix to zero + * and then keep the default value alpha=1. + * + * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply + * call this function with u.adjoint(). + */ + template<typename DerivedU> + SparseSelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1)); + + // const SparseLLT<PlainMatrixType, UpLo> llt() const; + // const SparseLDLT<PlainMatrixType, UpLo> ldlt() const; + + protected: + + const typename MatrixType::Nested m_matrix; +}; + +/*************************************************************************** +* Implementation of SparseMatrixBase methods +***************************************************************************/ + +template<typename Derived> +template<unsigned int UpLo> +const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const +{ + return derived(); +} + +template<typename Derived> +template<unsigned int UpLo> +SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() +{ + return derived(); +} + +/*************************************************************************** +* Implementation of SparseSelfAdjointView methods +***************************************************************************/ + +template<typename MatrixType, unsigned int UpLo> +template<typename DerivedU> +SparseSelfAdjointView<MatrixType,UpLo>& +SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha) +{ + SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint(); + if(alpha==Scalar(0)) + m_matrix = tmp.template triangularView<UpLo>(); + else + m_matrix += alpha * tmp.template triangularView<UpLo>(); + + return this; +} + +/*************************************************************************** +* Implementation of sparse self-adjoint time dense matrix +***************************************************************************/ + +template<typename Lhs, typename Rhs, int UpLo> +struct ei_traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> > + : ei_traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> > +{}; + +template<typename Lhs, typename Rhs, int UpLo> +class SparseSelfAdjointTimeDenseProduct + : public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct) + + SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const + { + // TODO use alpha + ei_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry"); + typedef typename ei_cleantype<Lhs>::type _Lhs; + typedef typename ei_cleantype<Rhs>::type _Rhs; + typedef typename _Lhs::InnerIterator LhsInnerIterator; + enum { + LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, + ProcessFirstHalf = + ((UpLo&(UpperTriangularBit|LowerTriangularBit))==(UpperTriangularBit|LowerTriangularBit)) + || ( (UpLo&UpperTriangularBit) && !LhsIsRowMajor) + || ( (UpLo&LowerTriangularBit) && LhsIsRowMajor), + ProcessSecondHalf = !ProcessFirstHalf + }; + for (int j=0; j<m_lhs.outerSize(); ++j) + { + LhsInnerIterator i(m_lhs,j); + if (ProcessSecondHalf && i && (i.index()==j)) + { + dest.row(j) += i.value() * m_rhs.row(j); + ++i; + } + Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0)); + for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) + { + int a = LhsIsRowMajor ? j : i.index(); + int b = LhsIsRowMajor ? i.index() : j; + typename Lhs::Scalar v = i.value(); + dest.row(a) += (v) * m_rhs.row(b); + dest.row(b) += ei_conj(v) * m_rhs.row(a); + } + if (ProcessFirstHalf && i && (i.index()==j)) + dest.row(j) += i.value() * m_rhs.row(j); + } + } + + private: + SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&); +}; + +template<typename Lhs, typename Rhs, int UpLo> +struct ei_traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> > + : ei_traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> > +{}; + +template<typename Lhs, typename Rhs, int UpLo> +class DenseTimeSparseSelfAdjointProduct + : public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct) + + DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const + { + // TODO + } + + private: + DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&); +}; +#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H diff --git a/Eigen/src/Sparse/SparseTriangular.h b/Eigen/src/Sparse/SparseTriangular.h deleted file mode 100644 index 42e7ff02a..000000000 --- a/Eigen/src/Sparse/SparseTriangular.h +++ /dev/null @@ -1,60 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SPARSE_TRIANGULAR_H -#define EIGEN_SPARSE_TRIANGULAR_H - -template<typename ExpressionType, int Mode> class SparseTriangular -{ - public: - - typedef typename ei_traits<ExpressionType>::Scalar Scalar; - typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret, - ExpressionType, const ExpressionType&>::ret ExpressionTypeNested; - typedef CwiseUnaryOp<ei_scalar_add_op<Scalar>, ExpressionType> ScalarAddReturnType; - - inline SparseTriangular(const ExpressionType& matrix) : m_matrix(matrix) {} - - /** \internal */ - inline const ExpressionType& _expression() const { return m_matrix; } - - template<typename OtherDerived> - typename ei_plain_matrix_type_column_major<OtherDerived>::type - solve(const MatrixBase<OtherDerived>& other) const; - template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const; - template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const; - - protected: - ExpressionTypeNested m_matrix; -}; - -template<typename Derived> -template<int Mode> -inline const SparseTriangular<Derived, Mode> -SparseMatrixBase<Derived>::triangular() const -{ - return derived(); -} - -#endif // EIGEN_SPARSE_TRIANGULAR_H diff --git a/Eigen/src/Sparse/SparseTriangularView.h b/Eigen/src/Sparse/SparseTriangularView.h new file mode 100644 index 000000000..6a9461528 --- /dev/null +++ b/Eigen/src/Sparse/SparseTriangularView.h @@ -0,0 +1,95 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_TRIANGULARVIEW_H +#define EIGEN_SPARSE_TRIANGULARVIEW_H + +template<typename MatrixType, int Mode> +struct ei_traits<SparseTriangularView<MatrixType,Mode> > +: public ei_traits<MatrixType> +{}; + +template<typename MatrixType, int Mode> class SparseTriangularView + : public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> > +{ + enum { SkipFirst = (Mode==LowerTriangular && (!MatrixType::Flags&RowMajorBit)) + || (Mode==UpperTriangular && ( MatrixType::Flags&RowMajorBit)) }; + public: + + class InnerIterator; + + inline int rows() { return m_matrix.rows(); } + inline int cols() { return m_matrix.cols(); } + + typedef typename ei_traits<MatrixType>::Scalar Scalar; + typedef typename ei_meta_if<ei_must_nest_by_value<MatrixType>::ret, + MatrixType, const MatrixType&>::ret MatrixTypeNested; + + inline SparseTriangularView(const MatrixType& matrix) : m_matrix(matrix) {} + + /** \internal */ + inline const MatrixType& nestedExpression() const { return m_matrix; } + + template<typename OtherDerived> + typename ei_plain_matrix_type_column_major<OtherDerived>::type + solve(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const; + template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const; + + protected: + MatrixTypeNested m_matrix; +}; + +template<typename MatrixType, int Mode> +class SparseTriangularView<MatrixType,Mode>::InnerIterator : public MatrixType::InnerIterator +{ + typedef typename MatrixType::InnerIterator Base; + public: + + EIGEN_STRONG_INLINE InnerIterator(const SparseTriangularView& view, int outer) + : Base(view.nestedExpression(), outer) + { + if(SkipFirst) + while((*this) && this->index()<outer) + ++(*this); + } + inline int row() const { return Base::row(); } + inline int col() const { return Base::col(); } + + EIGEN_STRONG_INLINE operator bool() const + { + return SkipFirst ? Base::operator bool() : (Base::operator bool() && this->index() < this->outer()); + } +}; + +template<typename Derived> +template<int Mode> +inline const SparseTriangularView<Derived, Mode> +SparseMatrixBase<Derived>::triangularView() const +{ + return derived(); +} + +#endif // EIGEN_SPARSE_TRIANGULARVIEW_H diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index b99be580c..6da1bee25 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -118,25 +118,29 @@ enum { }; template<typename Derived> class SparseMatrixBase; -template<typename _Scalar, int _Flags = 0> class SparseMatrix; -template<typename _Scalar, int _Flags = 0> class DynamicSparseMatrix; -template<typename _Scalar, int _Flags = 0> class SparseVector; -template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix; - -template<typename MatrixType> class SparseNestByValue; -template<typename MatrixType, int Size> class SparseInnerVectorSet; -template<typename ExpressionType, - unsigned int Added, unsigned int Removed> class SparseFlagged; -template<typename ExpressionType, int Mode> class SparseTriangular; -template<typename Lhs, typename Rhs> class SparseDiagonalProduct; +template<typename _Scalar, int _Flags = 0> class SparseMatrix; +template<typename _Scalar, int _Flags = 0> class DynamicSparseMatrix; +template<typename _Scalar, int _Flags = 0> class SparseVector; +template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix; + +template<typename MatrixType> class SparseNestByValue; +template<typename MatrixType, int Size> class SparseInnerVectorSet; +template<typename MatrixType, int Mode> class SparseTriangularView; +template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView; +template<typename Lhs, typename Rhs> class SparseDiagonalProduct; + + +template<typename Lhs, typename Rhs> class SparseProduct; +template<typename Lhs, typename Rhs> class SparseTimeDenseProduct; +template<typename Lhs, typename Rhs> class DenseTimeSparseProduct; template<typename Lhs, typename Rhs, typename LhsStorage = typename ei_traits<Lhs>::StorageType, typename RhsStorage = typename ei_traits<Rhs>::StorageType> struct ei_sparse_product_mode; -template<typename Lhs, typename Rhs, int ProductMode = ei_sparse_product_mode<Lhs,Rhs>::value> struct SparseProductReturnType; +template<typename Lhs, typename Rhs> struct SparseProductReturnType; -const int CoherentAccessPattern = 0x1; +const int CoherentAccessPattern = 0x1; const int InnerRandomAccessPattern = 0x2 | CoherentAccessPattern; const int OuterRandomAccessPattern = 0x4 | CoherentAccessPattern; const int RandomAccessPattern = 0x8 | OuterRandomAccessPattern | InnerRandomAccessPattern; diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h index cb2bc4d58..2c5b485b0 100644 --- a/Eigen/src/Sparse/TriangularSolver.h +++ b/Eigen/src/Sparse/TriangularSolver.h @@ -160,7 +160,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,ColMajor template<typename ExpressionType,int Mode> template<typename OtherDerived> -void SparseTriangular<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const +void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const { ei_assert(m_matrix.cols() == m_matrix.rows()); ei_assert(m_matrix.cols() == other.rows()); @@ -182,7 +182,7 @@ void SparseTriangular<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived template<typename ExpressionType,int Mode> template<typename OtherDerived> typename ei_plain_matrix_type_column_major<OtherDerived>::type -SparseTriangular<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const +SparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const { typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other); solveInPlace(res); @@ -210,10 +210,10 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> const bool IsLowerTriangular = (UpLo==LowerTriangular); AmbiVector<Scalar> tempVector(other.rows()*2); tempVector.setBounds(0,other.rows()); - + Rhs res(other.rows(), other.cols()); res.reserve(other.nonZeros()); - + for(int col=0 ; col<other.cols() ; ++col) { // FIXME estimate number of non zeros @@ -224,7 +224,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> { tempVector.coeffRef(rhsIt.index()) = rhsIt.value(); } - + for(int i=IsLowerTriangular?0:lhs.cols()-1; IsLowerTriangular?i<lhs.cols():i>=0; i+=IsLowerTriangular?1:-1) @@ -233,7 +233,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> Scalar& ci = tempVector.coeffRef(i); if (ci!=Scalar(0)) { - // find + // find typename Lhs::InnerIterator it(lhs, i); if(!(Mode & UnitDiagBit)) { @@ -260,8 +260,8 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> } } } - - + + int count = 0; // FIXME compute a reference value to filter zeros for (typename AmbiVector<Scalar>::Iterator it(tempVector/*,1e-12*/); it; ++it) @@ -281,7 +281,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> template<typename ExpressionType,int Mode> template<typename OtherDerived> -void SparseTriangular<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const +void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const { ei_assert(m_matrix.cols() == m_matrix.rows()); ei_assert(m_matrix.cols() == other.rows()); |