diff --git a/Eigen/Sparse b/Eigen/Sparse index 552dbde4a..f4dcad07e 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -98,7 +98,6 @@ struct Sparse {}; #include "src/Sparse/SparseVector.h" #include "src/Sparse/CoreIterators.h" #include "src/Sparse/SparseTranspose.h" -#include "src/Sparse/SparseCwise.h" #include "src/Sparse/SparseCwiseUnaryOp.h" #include "src/Sparse/SparseCwiseBinaryOp.h" #include "src/Sparse/SparseDot.h" @@ -107,12 +106,12 @@ struct Sparse {}; #include "src/Sparse/SparseFuzzy.h" #include "src/Sparse/SparseProduct.h" #include "src/Sparse/SparseDiagonalProduct.h" -#include "src/Sparse/SparseTriangular.h" +#include "src/Sparse/SparseTriangularView.h" +#include "src/Sparse/SparseSelfAdjointView.h" #include "src/Sparse/TriangularSolver.h" #include "src/Sparse/SparseLLT.h" #include "src/Sparse/SparseLDLT.h" #include "src/Sparse/SparseLU.h" -#include "src/Sparse/SparseExpressionMaker.h" #ifdef EIGEN_CHOLMOD_SUPPORT # include "src/Sparse/CholmodSupport.h" 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 class MatrixBase typedef Homogeneous::ColsAtCompileTime==1?Vertical:Horizontal> HomogeneousReturnType; const HomogeneousReturnType homogeneous() const; -/////////// Sparse module /////////// - - // dense = sparse * dense - template - Derived& lazyAssign(const SparseProduct& product); - // dense = dense * sparse - template - Derived& lazyAssign(const SparseProduct& 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 -template -const SelfAdjointView MatrixBase::selfadjointView() const +template +const SelfAdjointView MatrixBase::selfadjointView() const { return derived(); } template -template -SelfAdjointView MatrixBase::selfadjointView() +template +SelfAdjointView MatrixBase::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 class Translation; template class UniformScaling; template class Homogeneous; -// Sparse module: -template 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 -// struct ei_traits > -// { -// typedef typename ei_result_of< -// BinaryOp( -// typename Lhs::Scalar, -// typename Rhs::Scalar -// ) -// >::type Scalar; -// typedef typename ei_promote_storage_type::StorageType, -// typename ei_traits::StorageType>::ret StorageType; -// typedef typename Lhs::Nested LhsNested; -// typedef typename Rhs::Nested RhsNested; -// typedef typename ei_unref::type _LhsNested; -// typedef typename ei_unref::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::Cost -// }; -// }; - template<> struct ei_promote_storage_type { typedef Sparse ret; }; @@ -82,16 +53,9 @@ class CwiseBinaryOpImpl : public SparseMatrixBase > { public: - class InnerIterator; - typedef CwiseBinaryOp Derived; EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) -// typedef typename ei_traits::LhsNested LhsNested; -// typedef typename ei_traits::RhsNested RhsNested; -// typedef typename ei_unref::type _LhsNested; -// typedef typename ei_unref::type _RhsNested; - }; template::solveInPlace(MatrixBase &b) const return false; if (m_matrix.nonZeros()>0) // otherwise L==I - m_matrix.template triangular().solveInPlace(b); + m_matrix.template triangularView().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().solveInPlace(b); + m_matrix.transpose().template triangularView().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::solveInPlace(MatrixBase &b) const const int size = m_matrix.rows(); ei_assert(size==b.rows()); - m_matrix.template triangular().solveInPlace(b); + m_matrix.template triangularView().solveInPlace(b); // FIXME should be simply .adjoint() but it fails to compile... if (NumTraits::IsComplex) { CholMatrixType aux = m_matrix.conjugate(); - aux.transpose().template triangular().solveInPlace(b); + aux.transpose().template triangularView().solveInPlace(b); } else - m_matrix.transpose().template triangular().solveInPlace(b); + m_matrix.transpose().template triangularView().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::InnerIterator inline Scalar& valueRef() { return const_cast(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 class SparseMatrixBase : public AnyMatrixBase - inline Derived& operator=(const SparseProduct& product); + inline Derived& operator=(const SparseProduct& product); friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m) { @@ -337,12 +337,13 @@ template class SparseMatrixBase : public AnyMatrixBase friend - const typename SparseProductReturnType::Type + const DenseTimeSparseProduct operator*(const MatrixBase& lhs, const Derived& rhs) - { return typename SparseProductReturnType::Type(lhs.derived(),rhs); } + { return DenseTimeSparseProduct(lhs.derived(),rhs); } + // sparse * dense (returns a dense object) template - const typename SparseProductReturnType::Type + const SparseTimeDenseProduct operator*(const MatrixBase &other) const; template @@ -360,7 +361,10 @@ template class SparseMatrixBase : public AnyMatrixBase& other) const; template - inline const SparseTriangular triangular() const; + inline const SparseTriangularView triangularView() const; + + template inline const SparseSelfAdjointView selfadjointView() const; + template inline SparseSelfAdjointView selfadjointView(); template Scalar dot(const MatrixBase& other) const; template Scalar dot(const SparseMatrixBase& 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 struct ei_sparse_product_mode { enum { value = SparseTimeSparseProduct }; }; -template struct ei_sparse_product_mode { enum { value = DenseTimeSparseProduct }; }; -template struct ei_sparse_product_mode { enum { value = SparseTimeDenseProduct }; }; - -template -struct SparseProductReturnType -{ - typedef const typename ei_nested::type LhsNested; - typedef const typename ei_nested::type RhsNested; - - typedef SparseProduct Type; -}; - -// sparse product return type specialization template -struct SparseProductReturnType +struct SparseProductReturnType { typedef typename ei_traits::Scalar Scalar; enum { @@ -60,11 +46,11 @@ struct SparseProductReturnType SparseMatrix, const typename ei_nested::type>::ret RhsNested; - typedef SparseProduct Type; + typedef SparseProduct Type; }; -template -struct ei_traits > +template +struct ei_traits > { // clean the nested types: typedef typename ei_cleantype::type _LhsNested; @@ -85,7 +71,6 @@ struct ei_traits > MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime, EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit), - ResultIsSparse = ProductMode==SparseTimeSparseProduct, RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), @@ -96,16 +81,14 @@ struct ei_traits > CoeffReadCost = Dynamic }; - typedef typename ei_meta_if::ret StorageType; + typedef Sparse StorageType; - typedef typename ei_meta_if >, - MatrixBase > >::ret Base; + typedef SparseMatrixBase > Base; }; -template +template class SparseProduct : ei_no_assignment_operator, - public ei_traits >::Base + public ei_traits >::Base { public: @@ -256,38 +239,14 @@ struct ei_sparse_product_selector } }; -// NOTE eventually let's transpose one argument even in this case since it might be expensive if -// the result is not dense. -// template -// struct ei_sparse_product_selector -// { -// 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 -// template -// inline Derived& SparseMatrixBase::lazyAssign(const SparseProduct& product) -// { -// // std::cout << "sparse product to dense\n"; -// ei_sparse_product_selector< -// typename ei_cleantype::type, -// typename ei_cleantype::type, -// typename ei_cleantype::type>::run(product.lhs(),product.rhs(),derived()); -// return derived(); -// } - // sparse = sparse * sparse template template -inline Derived& SparseMatrixBase::operator=(const SparseProduct& product) +inline Derived& SparseMatrixBase::operator=(const SparseProduct& product) { ei_sparse_product_selector< typename ei_cleantype::type, @@ -296,78 +255,79 @@ inline Derived& SparseMatrixBase::operator=(const SparseProduct -EIGEN_DONT_INLINE void ei_sparse_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) -{ - typedef typename ei_cleantype::type _Lhs; - typedef typename ei_cleantype::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 dst_j(dst.row(LhsIsRowMajor ? j : 0)); - for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) - { - if(LhsIsSelfAdjoint) - { - 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); - } - 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 template -Derived& MatrixBase::lazyAssign(const SparseProduct& product) +struct ei_traits > + : ei_traits, Lhs, Rhs> > { - derived().setZero(); - ei_sparse_time_dense_product(product.lhs(), product.rhs(), derived()); - return derived(); -} + typedef Dense StorageType; +}; + +template +class SparseTimeDenseProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseTimeDenseProduct) + + SparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template void scaleAndAddTo(Dest& dest, Scalar alpha) const + { + typedef typename ei_cleantype::type _Lhs; + typedef typename ei_cleantype::type _Rhs; + typedef typename _Lhs::InnerIterator LhsInnerIterator; + enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit }; + for(int j=0; j 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); + } + } + } + + private: + SparseTimeDenseProduct& operator=(const SparseTimeDenseProduct&); +}; + // dense = dense * sparse -template template -Derived& MatrixBase::lazyAssign(const SparseProduct& product) +struct ei_traits > + : ei_traits, Lhs, Rhs> > { - typedef typename ei_cleantype::type _Rhs; - typedef typename _Rhs::InnerIterator RhsInnerIterator; - enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit }; - derived().setZero(); - for (int j=0; j +class DenseTimeSparseProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseProduct) + + DenseTimeSparseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template void scaleAndAddTo(Dest& dest, Scalar alpha) const + { + typedef typename ei_cleantype::type _Rhs; + typedef typename _Rhs::InnerIterator RhsInnerIterator; + enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit }; + for(int j=0; j @@ -381,10 +341,10 @@ SparseMatrixBase::operator*(const SparseMatrixBase &other // sparse * dense template template -EIGEN_STRONG_INLINE const typename SparseProductReturnType::Type +EIGEN_STRONG_INLINE const SparseTimeDenseProduct SparseMatrixBase::operator*(const MatrixBase &other) const { - return typename SparseProductReturnType::Type(derived(), other.derived()); + return SparseTimeDenseProduct(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 +// +// 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 . + +#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 +class SparseSelfAdjointTimeDenseProduct; + +template +class DenseTimeSparseSelfAdjointProduct; + +template class SparseSelfAdjointView +{ + public: + + typedef typename ei_traits::Scalar Scalar; + + inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix) + { + ei_assert(ei_are_flags_consistent::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 + SparseSelfAdjointTimeDenseProduct + operator*(const MatrixBase& rhs) const + { + return SparseSelfAdjointTimeDenseProduct(m_matrix, rhs.derived()); + } + + /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ + template friend + DenseTimeSparseSelfAdjointProduct + operator*(const MatrixBase& lhs, const SparseSelfAdjointView& rhs) + { + return DenseTimeSparseSelfAdjointProduct(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 + SparseSelfAdjointView& rankUpdate(const MatrixBase& u, Scalar alpha = Scalar(1)); + + // const SparseLLT llt() const; + // const SparseLDLT ldlt() const; + + protected: + + const typename MatrixType::Nested m_matrix; +}; + +/*************************************************************************** +* Implementation of SparseMatrixBase methods +***************************************************************************/ + +template +template +const SparseSelfAdjointView SparseMatrixBase::selfadjointView() const +{ + return derived(); +} + +template +template +SparseSelfAdjointView SparseMatrixBase::selfadjointView() +{ + return derived(); +} + +/*************************************************************************** +* Implementation of SparseSelfAdjointView methods +***************************************************************************/ + +template +template +SparseSelfAdjointView& +SparseSelfAdjointView::rankUpdate(const MatrixBase& u, Scalar alpha) +{ + SparseMatrix tmp = u * u.adjoint(); + if(alpha==Scalar(0)) + m_matrix = tmp.template triangularView(); + else + m_matrix += alpha * tmp.template triangularView(); + + return this; +} + +/*************************************************************************** +* Implementation of sparse self-adjoint time dense matrix +***************************************************************************/ + +template +struct ei_traits > + : ei_traits, Lhs, Rhs> > +{}; + +template +class SparseSelfAdjointTimeDenseProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct) + + SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template 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::type _Lhs; + typedef typename ei_cleantype::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 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 +struct ei_traits > + : ei_traits, Lhs, Rhs> > +{}; + +template +class DenseTimeSparseSelfAdjointProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct) + + DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template 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 -// -// 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 . - -#ifndef EIGEN_SPARSE_TRIANGULAR_H -#define EIGEN_SPARSE_TRIANGULAR_H - -template class SparseTriangular -{ - public: - - typedef typename ei_traits::Scalar Scalar; - typedef typename ei_meta_if::ret, - ExpressionType, const ExpressionType&>::ret ExpressionTypeNested; - typedef CwiseUnaryOp, ExpressionType> ScalarAddReturnType; - - inline SparseTriangular(const ExpressionType& matrix) : m_matrix(matrix) {} - - /** \internal */ - inline const ExpressionType& _expression() const { return m_matrix; } - - template - typename ei_plain_matrix_type_column_major::type - solve(const MatrixBase& other) const; - template void solveInPlace(MatrixBase& other) const; - template void solveInPlace(SparseMatrixBase& other) const; - - protected: - ExpressionTypeNested m_matrix; -}; - -template -template -inline const SparseTriangular -SparseMatrixBase::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 +// +// 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 . + +#ifndef EIGEN_SPARSE_TRIANGULARVIEW_H +#define EIGEN_SPARSE_TRIANGULARVIEW_H + +template +struct ei_traits > +: public ei_traits +{}; + +template class SparseTriangularView + : public SparseMatrixBase > +{ + 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::Scalar Scalar; + typedef typename ei_meta_if::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 ei_plain_matrix_type_column_major::type + solve(const MatrixBase& other) const; + + template void solveInPlace(MatrixBase& other) const; + template void solveInPlace(SparseMatrixBase& other) const; + + protected: + MatrixTypeNested m_matrix; +}; + +template +class SparseTriangularView::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()index() < this->outer()); + } +}; + +template +template +inline const SparseTriangularView +SparseMatrixBase::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 class SparseMatrixBase; -template class SparseMatrix; -template class DynamicSparseMatrix; -template class SparseVector; -template class MappedSparseMatrix; +template class SparseMatrix; +template class DynamicSparseMatrix; +template class SparseVector; +template class MappedSparseMatrix; -template class SparseNestByValue; -template class SparseInnerVectorSet; -template class SparseFlagged; -template class SparseTriangular; -template class SparseDiagonalProduct; +template class SparseNestByValue; +template class SparseInnerVectorSet; +template class SparseTriangularView; +template class SparseSelfAdjointView; +template class SparseDiagonalProduct; + + +template class SparseProduct; +template class SparseTimeDenseProduct; +template class DenseTimeSparseProduct; template::StorageType, typename RhsStorage = typename ei_traits::StorageType> struct ei_sparse_product_mode; -template::value> struct SparseProductReturnType; +template 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 template -void SparseTriangular::solveInPlace(MatrixBase& other) const +void SparseTriangularView::solveInPlace(MatrixBase& other) const { ei_assert(m_matrix.cols() == m_matrix.rows()); ei_assert(m_matrix.cols() == other.rows()); @@ -182,7 +182,7 @@ void SparseTriangular::solveInPlace(MatrixBase template typename ei_plain_matrix_type_column_major::type -SparseTriangular::solve(const MatrixBase& other) const +SparseTriangularView::solve(const MatrixBase& other) const { typename ei_plain_matrix_type_column_major::type res(other); solveInPlace(res); @@ -210,10 +210,10 @@ struct ei_sparse_solve_triangular_sparse_selector const bool IsLowerTriangular = (UpLo==LowerTriangular); AmbiVector tempVector(other.rows()*2); tempVector.setBounds(0,other.rows()); - + Rhs res(other.rows(), other.cols()); res.reserve(other.nonZeros()); - + for(int col=0 ; col { tempVector.coeffRef(rhsIt.index()) = rhsIt.value(); } - + for(int i=IsLowerTriangular?0:lhs.cols()-1; IsLowerTriangular?i=0; i+=IsLowerTriangular?1:-1) @@ -233,7 +233,7 @@ struct ei_sparse_solve_triangular_sparse_selector 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 } } } - - + + int count = 0; // FIXME compute a reference value to filter zeros for (typename AmbiVector::Iterator it(tempVector/*,1e-12*/); it; ++it) @@ -281,7 +281,7 @@ struct ei_sparse_solve_triangular_sparse_selector template template -void SparseTriangular::solveInPlace(SparseMatrixBase& other) const +void SparseTriangularView::solveInPlace(SparseMatrixBase& other) const { ei_assert(m_matrix.cols() == m_matrix.rows()); ei_assert(m_matrix.cols() == other.rows()); diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 3f0e793d5..e944d6c53 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -114,10 +114,10 @@ template void sparse_product(const SparseMatrixType& VERIFY_IS_APPROX(mS.transpose().conjugate(), mS); VERIFY_IS_APPROX(mS, refS); VERIFY_IS_APPROX(x=mS*b, refX=refS*b); - // TODO properly implement triangular/selfadjoint views -// VERIFY_IS_APPROX(x=mUp.template marked()*b, refX=refS*b); -// VERIFY_IS_APPROX(x=mLo.template marked()*b, refX=refS*b); -// VERIFY_IS_APPROX(x=mS.template marked()*b, refX=refS*b); + + VERIFY_IS_APPROX(x=mUp.template selfadjointView()*b, refX=refS*b); + VERIFY_IS_APPROX(x=mLo.template selfadjointView()*b, refX=refS*b); + VERIFY_IS_APPROX(x=mS.template selfadjointView()*b, refX=refS*b); } } diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp index a530a9515..24107977c 100644 --- a/test/sparse_solvers.cpp +++ b/test/sparse_solvers.cpp @@ -66,12 +66,12 @@ template void sparse_solvers(int rows, int cols) // lower - dense initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.template triangularView().solve(vec2), - m2.template triangular().solve(vec3)); + m2.template triangularView().solve(vec3)); // upper - dense initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.template triangularView().solve(vec2), - m2.template triangular().solve(vec3)); + m2.template triangularView().solve(vec3)); // TODO test row major @@ -82,20 +82,20 @@ template void sparse_solvers(int rows, int cols) initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular); initSparse(density, refMatB, matB); refMat2.template triangularView().solveInPlace(refMatB); - m2.template triangular().solveInPlace(matB); + m2.template triangularView().solveInPlace(matB); VERIFY_IS_APPROX(matB.toDense(), refMatB); // upper - sparse initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular); initSparse(density, refMatB, matB); refMat2.template triangularView().solveInPlace(refMatB); - m2.template triangular().solveInPlace(matB); + m2.template triangularView().solveInPlace(matB); VERIFY_IS_APPROX(matB, refMatB); // test deprecated API initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.template triangularView().solve(vec2), - m2.template triangular().solve(vec3)); + m2.template triangularView().solve(vec3)); } // test LLT