diff options
| -rw-r--r-- | Eigen/Sparse | 1 | ||||
| -rw-r--r-- | Eigen/src/Sparse/AmbiVector.h | 15 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseBlock.h | 36 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseLDLT.h | 8 | ||||
| -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 | 7 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseTriangular.h | 60 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SparseUtil.h | 1 | ||||
| -rw-r--r-- | Eigen/src/Sparse/SuperLUSupport.h | 2 | ||||
| -rw-r--r-- | Eigen/src/Sparse/TriangularSolver.h | 218 | ||||
| -rw-r--r-- | test/sparse_product.cpp | 6 | ||||
| -rw-r--r-- | test/sparse_solvers.cpp | 32 |
13 files changed, 332 insertions, 61 deletions
diff --git a/Eigen/Sparse b/Eigen/Sparse index 6de24d34d..b0f2b2a12 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -104,6 +104,7 @@ namespace Eigen { #include "src/Sparse/SparseFlagged.h" #include "src/Sparse/SparseProduct.h" #include "src/Sparse/SparseDiagonalProduct.h" +#include "src/Sparse/SparseTriangular.h" #include "src/Sparse/TriangularSolver.h" #include "src/Sparse/SparseLLT.h" #include "src/Sparse/SparseLDLT.h" diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h index 75001a2fa..e7a7c8784 100644 --- a/Eigen/src/Sparse/AmbiVector.h +++ b/Eigen/src/Sparse/AmbiVector.h @@ -36,7 +36,7 @@ template<typename _Scalar> class AmbiVector typedef _Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; AmbiVector(int size) - : m_buffer(0), m_size(0), m_allocatedSize(0), m_allocatedElements(0), m_mode(-1) + : m_buffer(0), m_zero(0), m_size(0), m_allocatedSize(0), m_allocatedElements(0), m_mode(-1) { resize(size); } @@ -44,7 +44,7 @@ template<typename _Scalar> class AmbiVector void init(RealScalar estimatedDensity); void init(int mode); - void nonZeros() const; + int nonZeros() const; /** Specifies a sub-vector to work on */ void setBounds(int start, int end) { m_start = start; m_end = end; } @@ -53,7 +53,7 @@ template<typename _Scalar> class AmbiVector void restart(); Scalar& coeffRef(int i); - Scalar coeff(int i); + Scalar& coeff(int i); class Iterator; @@ -112,6 +112,7 @@ template<typename _Scalar> class AmbiVector // used to store data in both mode Scalar* m_buffer; + Scalar m_zero; int m_size; int m_start; int m_end; @@ -131,7 +132,7 @@ template<typename _Scalar> class AmbiVector /** \returns the number of non zeros in the current sub vector */ template<typename Scalar> -void AmbiVector<Scalar>::nonZeros() const +int AmbiVector<Scalar>::nonZeros() const { if (m_mode==IsSparse) return m_llSize; @@ -254,7 +255,7 @@ Scalar& AmbiVector<Scalar>::coeffRef(int i) } template<typename Scalar> -Scalar AmbiVector<Scalar>::coeff(int i) +Scalar& AmbiVector<Scalar>::coeff(int i) { if (m_mode==IsDense) return m_buffer[i]; @@ -264,7 +265,7 @@ Scalar AmbiVector<Scalar>::coeff(int i) ei_assert(m_mode==IsSparse); if ((m_llSize==0) || (i<llElements[m_llStart].index)) { - return Scalar(0); + return m_zero; } else { @@ -275,7 +276,7 @@ Scalar AmbiVector<Scalar>::coeff(int i) if (llElements[elid].index==i) return llElements[m_llCurrent].value; else - return Scalar(0); + return m_zero; } } } diff --git a/Eigen/src/Sparse/SparseBlock.h b/Eigen/src/Sparse/SparseBlock.h index c39066676..86567f1a5 100644 --- a/Eigen/src/Sparse/SparseBlock.h +++ b/Eigen/src/Sparse/SparseBlock.h @@ -150,6 +150,21 @@ class SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size> { return operator=<SparseInnerVectorSet>(other); } + + int nonZeros() const + { + int count = 0; + for (int j=0; j<m_outerSize; ++j) + count += m_matrix._data()[m_outerStart+j].size(); + return count; + } + + const Scalar& lastCoeff() const + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(SparseInnerVectorSet); + ei_assert(m_matrix.data()[m_outerStart].size()>0); + return m_matrix.data()[m_outerStart].vale(m_matrix.data()[m_outerStart].size()-1); + } // template<typename Sparse> // inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other) @@ -172,12 +187,12 @@ class SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size> /*************************************************************************** * specialisation for SparseMatrix ***************************************************************************/ -/* + template<typename _Scalar, int _Options, int Size> class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> : public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> > { - typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType; + typedef SparseMatrix<_Scalar, _Options> MatrixType; enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor }; public: @@ -233,7 +248,20 @@ class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> { return m_matrix._valuePtr() + m_matrix._outerIndexPtr()[m_outerStart]; } inline const int* _innerIndexPtr() const { return m_matrix._innerIndexPtr() + m_matrix._outerIndexPtr()[m_outerStart]; } - inline const int* _outerIndexPtr() const { return m_matrix._outerIndexPtr() + m_outerStart; } + inline const int* _outerIndexPtr() const + { return m_matrix._outerIndexPtr() + m_outerStart; } + + int nonZeros() const + { + return size_t(m_matrix._outerIndexPtr()[m_outerStart+m_outerSize.value()]) + - size_t(m_matrix._outerIndexPtr()[m_outerStart]); } + + const Scalar& lastCoeff() const + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(SparseInnerVectorSet); + ei_assert(nonZeros()>0); + return m_matrix._valuePtr()[m_matrix._outerIndexPtr()[m_outerStart+1]-1]; + } // template<typename Sparse> // inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other) @@ -251,7 +279,7 @@ class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> const ei_int_if_dynamic<Size> m_outerSize; }; -*/ + //---------- /** \returns the i-th row of the matrix \c *this. For row-major matrix only. */ diff --git a/Eigen/src/Sparse/SparseLDLT.h b/Eigen/src/Sparse/SparseLDLT.h index a1bac4d08..80d67ec8a 100644 --- a/Eigen/src/Sparse/SparseLDLT.h +++ b/Eigen/src/Sparse/SparseLDLT.h @@ -79,7 +79,7 @@ class SparseLDLT protected: typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - typedef SparseMatrix<Scalar,LowerTriangular|UnitDiagBit> CholMatrixType; + typedef SparseMatrix<Scalar> CholMatrixType; typedef Matrix<Scalar,MatrixType::ColsAtCompileTime,1> VectorType; enum { @@ -137,7 +137,7 @@ class SparseLDLT * overloads the MemoryEfficient flags) * * \sa flags() */ - void settagss(int f) { m_flags = f; } + void settags(int f) { m_flags = f; } /** \returns the current flags */ int flags() const { return m_flags; } @@ -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.solveTriangularInPlace(b); + m_matrix.template triangular<LowerTriangular|UnitDiagBit>().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().solveTriangularInPlace(b); + m_matrix.transpose().template triangular<UpperTriangular|UnitDiagBit>().solveInPlace(b); return true; } diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h index 864c4415a..40cefa189 100644 --- a/Eigen/src/Sparse/SparseLLT.h +++ b/Eigen/src/Sparse/SparseLLT.h @@ -191,15 +191,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 triangular<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 triangular<UpperTriangular>().solveInPlace(b); } else - m_matrix.transpose().template triangular<UpperTriangular>.solveInPlace(b); + m_matrix.transpose().template triangular<UpperTriangular>().solveInPlace(b); return true; } diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 3f09596bc..ad94e63ba 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -164,6 +164,7 @@ class SparseMatrix { ei_assert(m_data.index(m_data.size()-1)<inner && "wrong sorted insertion"); } +// std::cerr << size_t(m_outerIndex[outer+1]) << " == " << m_data.size() << "\n"; assert(size_t(m_outerIndex[outer+1]) == m_data.size()); int id = m_outerIndex[outer+1]; ++m_outerIndex[outer+1]; diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index f20841ce1..e611744d4 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -308,12 +308,19 @@ template<typename Derived> class SparseMatrixBase template<typename OtherDerived> Derived& operator*=(const SparseMatrixBase<OtherDerived>& other); + // deprecated template<typename OtherDerived> typename ei_plain_matrix_type_column_major<OtherDerived>::type solveTriangular(const MatrixBase<OtherDerived>& other) const; + // deprecated template<typename OtherDerived> void solveTriangularInPlace(MatrixBase<OtherDerived>& other) const; +// template<typename OtherDerived> +// void solveTriangularInPlace(SparseMatrixBase<OtherDerived>& other) const; + + template<int Mode> + inline const SparseTriangular<Derived, Mode> triangular() const; 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/SparseTriangular.h b/Eigen/src/Sparse/SparseTriangular.h new file mode 100644 index 000000000..ca16caf96 --- /dev/null +++ b/Eigen/src/Sparse/SparseTriangular.h @@ -0,0 +1,60 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// 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/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index 393cdda6e..6610bdf30 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -113,6 +113,7 @@ template<typename UnaryOp, typename MatrixType> class SparseCwiseUnaryO template<typename BinaryOp, typename Lhs, typename Rhs> class SparseCwiseBinaryOp; 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 Lhs, typename Rhs> struct ei_sparse_product_mode; diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h index 3c9a4fcce..cf0f802ed 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/Eigen/src/Sparse/SuperLUSupport.h @@ -543,7 +543,7 @@ typename SparseLU<MatrixType,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::dete if (m_extractedDataAreDirty) extractData(); - // TODO this code coule be moved to the default/base backend + // TODO this code could be moved to the default/base backend // FIXME perhaps we have to take into account the scale factors m_sluRscale and m_sluCscale ??? Scalar det = Scalar(1); for (int j=0; j<m_u.cols(); ++j) diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h index 8948ae45e..168439b30 100644 --- a/Eigen/src/Sparse/TriangularSolver.h +++ b/Eigen/src/Sparse/TriangularSolver.h @@ -25,9 +25,18 @@ #ifndef EIGEN_SPARSETRIANGULARSOLVER_H #define EIGEN_SPARSETRIANGULARSOLVER_H +template<typename Lhs, typename Rhs, int Mode, + int UpLo = (Mode & LowerTriangularBit) + ? LowerTriangular + : (Mode & UpperTriangularBit) + ? UpperTriangular + : -1, + int StorageOrder = int(ei_traits<Lhs>::Flags) & RowMajorBit> +struct ei_sparse_solve_triangular_selector; + // forward substitution, row-major -template<typename Lhs, typename Rhs> -struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,RowMajor|IsSparse> +template<typename Lhs, typename Rhs, int Mode> +struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,RowMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -45,7 +54,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,RowMajor|IsSparse> lastIndex = it.index(); tmp -= lastVal * other.coeff(lastIndex,col); } - if (Lhs::Flags & UnitDiagBit) + if (Mode & UnitDiagBit) other.coeffRef(i,col) = tmp; else { @@ -58,8 +67,8 @@ struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,RowMajor|IsSparse> }; // backward substitution, row-major -template<typename Lhs, typename Rhs> -struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,RowMajor|IsSparse> +template<typename Lhs, typename Rhs, int Mode> +struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,RowMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -77,7 +86,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,RowMajor|IsSparse> tmp -= it.value() * other.coeff(it.index(),col); } - if (Lhs::Flags & UnitDiagBit) + if (Mode & UnitDiagBit) other.coeffRef(i,col) = tmp; else { @@ -91,8 +100,8 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,RowMajor|IsSparse> }; // forward substitution, col-major -template<typename Lhs, typename Rhs> -struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,ColMajor|IsSparse> +template<typename Lhs, typename Rhs, int Mode> +struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,LowerTriangular,ColMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -101,26 +110,28 @@ struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,ColMajor|IsSparse> { for(int i=0; i<lhs.cols(); ++i) { - typename Lhs::InnerIterator it(lhs, i); - if(!(Lhs::Flags & UnitDiagBit)) + Scalar& tmp = other.coeffRef(i,col); + if (tmp!=Scalar(0)) // optimization when other is actually sparse { - // std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n"; - ei_assert(it.index()==i); - other.coeffRef(i,col) /= it.value(); + typename Lhs::InnerIterator it(lhs, i); + if(!(Mode & UnitDiagBit)) + { + ei_assert(it.index()==i); + tmp /= it.value(); + } + if (it.index()==i) + ++it; + for(; it; ++it) + other.coeffRef(it.index(), col) -= tmp * it.value(); } - Scalar tmp = other.coeffRef(i,col); - if (it.index()==i) - ++it; - for(; it; ++it) - other.coeffRef(it.index(), col) -= tmp * it.value(); } } } }; // backward substitution, col-major -template<typename Lhs, typename Rhs> -struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,ColMajor|IsSparse> +template<typename Lhs, typename Rhs, int Mode> +struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,ColMajor> { typedef typename Rhs::Scalar Scalar; static void run(const Lhs& lhs, Rhs& other) @@ -129,29 +140,32 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,ColMajor|IsSparse> { for(int i=lhs.cols()-1; i>=0; --i) { - if(!(Lhs::Flags & UnitDiagBit)) + Scalar& tmp = other.coeffRef(i,col); + if (tmp!=Scalar(0)) // optimization when other is actually sparse { - // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the - // last element of the column ! - other.coeffRef(i,col) /= lhs.coeff(i,i); + if(!(Mode & UnitDiagBit)) + { + // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the + // last element of the column ! + other.coeffRef(i,col) /= lhs.innerVector(i).lastCoeff(); + } + typename Lhs::InnerIterator it(lhs, i); + for(; it && it.index()<i; ++it) + other.coeffRef(it.index(), col) -= tmp * it.value(); } - Scalar tmp = other.coeffRef(i,col); - typename Lhs::InnerIterator it(lhs, i); - for(; it && it.index()<i; ++it) - other.coeffRef(it.index(), col) -= tmp * it.value(); } } } }; -template<typename Derived> +template<typename ExpressionType,int Mode> template<typename OtherDerived> -void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const +void SparseTriangular<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const { - ei_assert(derived().cols() == derived().rows()); - ei_assert(derived().cols() == other.rows()); - ei_assert(!(Flags & ZeroDiagBit)); - ei_assert(Flags & (UpperTriangularBit|LowerTriangularBit)); + ei_assert(m_matrix.cols() == m_matrix.rows()); + ei_assert(m_matrix.cols() == other.rows()); + ei_assert(!(Mode & ZeroDiagBit)); + ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit }; @@ -159,19 +173,151 @@ void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& typename ei_plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy; OtherCopy otherCopy(other.derived()); - ei_solve_triangular_selector<Derived, typename ei_unref<OtherCopy>::type>::run(derived(), otherCopy); + ei_sparse_solve_triangular_selector<ExpressionType, typename ei_unref<OtherCopy>::type, Mode>::run(m_matrix, otherCopy); if (copy) other = otherCopy; } +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 +{ + typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other); + solveInPlace(res); + return res; +} + +// pure sparse path + +template<typename Lhs, typename Rhs, int Mode, + int UpLo = (Mode & LowerTriangularBit) + ? LowerTriangular + : (Mode & UpperTriangularBit) + ? UpperTriangular + : -1, + int StorageOrder = int(Lhs::Flags) & (RowMajorBit)> +struct ei_sparse_solve_triangular_sparse_selector; + +// forward substitution, col-major +template<typename Lhs, typename Rhs, int Mode, int UpLo> +struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + const bool IsLowerTriangular = (UpLo==LowerTriangular); + AmbiVector<Scalar> tempVector(other.rows()*2); + tempVector.setBounds(0,other.rows()); + + Rhs res(other.rows(), other.cols()); + res.startFill(other.nonZeros()); + + for(int col=0 ; col<other.cols() ; ++col) + { + // FIXME estimate number of non zeros + tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/); + tempVector.setZero(); + tempVector.restart(); + for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt) + { + tempVector.coeffRef(rhsIt.index()) = rhsIt.value(); + } + + for(int i=IsLowerTriangular?0:lhs.cols()-1; + IsLowerTriangular?i<lhs.cols():i>=0; + i+=IsLowerTriangular?1:-1) + { + tempVector.restart(); + Scalar& ci = tempVector.coeffRef(i); + if (ci!=Scalar(0)) + { + // find + typename Lhs::InnerIterator it(lhs, i); + if(!(Mode & UnitDiagBit)) + { + if (IsLowerTriangular) + { + ei_assert(it.index()==i); + ci /= it.value(); + } + else + ci /= lhs.coeff(i,i); + } + tempVector.restart(); + if (IsLowerTriangular) + { + if (it.index()==i) + ++it; + for(; it; ++it) + tempVector.coeffRef(it.index()) -= ci * it.value(); + } + else + { + for(; it && it.index()<i; ++it) + tempVector.coeffRef(it.index()) -= ci * it.value(); + } + } + } + + + int count = 0; + // FIXME compute a reference value to filter zeros + for (typename AmbiVector<Scalar>::Iterator it(tempVector/*,1e-12*/); it; ++it) + { + ++ count; +// std::cerr << "fill " << it.index() << ", " << col << "\n"; +// std::cout << it.value() << " "; + res.fill(it.index(), col) = it.value(); + } +// std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n"; + } + res.endFill(); + other = res.markAsRValue(); + } +}; + +template<typename ExpressionType,int Mode> +template<typename OtherDerived> +void SparseTriangular<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const +{ + ei_assert(m_matrix.cols() == m_matrix.rows()); + ei_assert(m_matrix.cols() == other.rows()); + ei_assert(!(Mode & ZeroDiagBit)); + ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit)); + +// enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit }; + +// typedef typename ei_meta_if<copy, +// typename ei_plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy; +// OtherCopy otherCopy(other.derived()); + + ei_sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(m_matrix, other.derived()); + +// if (copy) +// other = otherCopy; +} + + +// deprecated stuff: + +/** \deprecated */ +template<typename Derived> +template<typename OtherDerived> +void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const +{ + this->template triangular<Flags&(UpperTriangularBit|LowerTriangularBit)>().solveInPlace(other); +} + +/** \deprecated */ template<typename Derived> template<typename OtherDerived> typename ei_plain_matrix_type_column_major<OtherDerived>::type SparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const { typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other); - solveTriangularInPlace(res); + derived().solveTriangularInPlace(res); return res; } diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 3a0356f00..00b5a3ed5 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -40,6 +40,7 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows); DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows); + DenseMatrix refMat5 = DenseMatrix::Random(rows, rows); DenseMatrix dm4 = DenseMatrix::Zero(rows, rows); SparseMatrixType m2(rows, rows); SparseMatrixType m3(rows, rows); @@ -57,7 +58,10 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose()); VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3); VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); - + + VERIFY_IS_APPROX(dm4=m2*(refMat3+refMat3), refMat4=refMat2*(refMat3+refMat3)); + VERIFY_IS_APPROX(dm4=m2.transpose()*(refMat3+refMat5)*0.5, refMat4=refMat2.transpose()*(refMat3+refMat5)*0.5); + // dense * sparse VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3); VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp index 8ce4e2264..d1090dfed 100644 --- a/test/sparse_solvers.cpp +++ b/test/sparse_solvers.cpp @@ -63,17 +63,39 @@ template<typename Scalar> void sparse_solvers(int rows, int cols) SparseMatrix<Scalar> m2(rows, cols); DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); - // lower + // lower - dense initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2), - m2.template marked<LowerTriangular>().solveTriangular(vec3)); + m2.template triangular<LowerTriangular>().solve(vec3)); - // upper + // upper - dense initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().solveTriangular(vec2), - m2.template marked<UpperTriangular>().solveTriangular(vec3)); - + m2.template triangular<UpperTriangular>().solve(vec3)); + // TODO test row major + + SparseMatrix<Scalar> matB(rows, rows); + DenseMatrix refMatB = DenseMatrix::Zero(rows, rows); + + // lower - sparse + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular); + initSparse<Scalar>(density, refMatB, matB); + refMat2.template marked<LowerTriangular>().solveTriangularInPlace(refMatB); + m2.template triangular<LowerTriangular>().solveInPlace(matB); + VERIFY_IS_APPROX(matB.toDense(), refMatB); + + // upper - sparse + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular); + initSparse<Scalar>(density, refMatB, matB); + refMat2.template marked<UpperTriangular>().solveTriangularInPlace(refMatB); + m2.template triangular<UpperTriangular>().solveInPlace(matB); + VERIFY_IS_APPROX(matB, refMatB); + + // test deprecated API + initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); + VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2), + m2.template marked<LowerTriangular>().solveTriangular(vec3)); } // test LLT |