Re-enable products with triangular views of sparse matrices: we simply have to treat them as a sparse matrix.
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 36bbd46..263d544 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h
@@ -247,7 +247,7 @@ inline const AdjointReturnType adjoint() const { return AdjointReturnType(m_matrix.adjoint()); } - typedef TriangularView<Transpose<MatrixType>,TransposeMode> TransposeReturnType; + typedef TriangularView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType; /** \sa MatrixBase::transpose() */ EIGEN_DEVICE_FUNC inline TransposeReturnType transpose() @@ -255,11 +255,13 @@ EIGEN_STATIC_ASSERT_LVALUE(MatrixType) return TransposeReturnType(m_matrix.const_cast_derived().transpose()); } + + typedef TriangularView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType; /** \sa MatrixBase::transpose() const */ EIGEN_DEVICE_FUNC - inline const TransposeReturnType transpose() const + inline const ConstTransposeReturnType transpose() const { - return TransposeReturnType(m_matrix.transpose()); + return ConstTransposeReturnType(m_matrix.transpose()); } template<typename Other>
diff --git a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h index 19d9eaa..a30522f 100644 --- a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h +++ b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h
@@ -141,6 +141,8 @@ typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrixAux; typedef typename sparse_eval<ColMajorMatrixAux,ResultType::RowsAtCompileTime,ResultType::ColsAtCompileTime>::type ColMajorMatrix; + // If the result is tall and thin (in the extreme case a column vector) + // then it is faster to sort the coefficients inplace instead of transposing twice. // FIXME, the following heuristic is probably not very good. if(lhs.rows()>=rhs.cols()) {
diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h index ed3cb19..5aea114 100644 --- a/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/Eigen/src/SparseCore/SparseDenseProduct.h
@@ -147,6 +147,11 @@ }; template<typename Lhs, typename Rhs, int ProductType> +struct generic_product_impl<Lhs, Rhs, SparseTriangularShape, DenseShape, ProductType> + : generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType> +{}; + +template<typename Lhs, typename Rhs, int ProductType> struct generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType> { template<typename Dest> @@ -158,12 +163,17 @@ RhsNested rhsNested(rhs); dst.setZero(); - // transpoe everything + // transpose everything Transpose<Dest> dstT(dst); internal::sparse_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, typename Dest::Scalar(1)); } }; +template<typename Lhs, typename Rhs, int ProductType> +struct generic_product_impl<Lhs, Rhs, DenseShape, SparseTriangularShape, ProductType> + : generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType> +{}; + template<typename LhsT, typename RhsT, bool NeedToTranspose> struct sparse_dense_outer_product_evaluator {
diff --git a/Eigen/src/SparseCore/SparseProduct.h b/Eigen/src/SparseCore/SparseProduct.h index ae70737..c62386e 100644 --- a/Eigen/src/SparseCore/SparseProduct.h +++ b/Eigen/src/SparseCore/SparseProduct.h
@@ -33,6 +33,7 @@ namespace internal { +// sparse * sparse template<typename Lhs, typename Rhs, int ProductType> struct generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType> { @@ -48,6 +49,18 @@ } }; +// sparse * sparse-triangular +template<typename Lhs, typename Rhs, int ProductType> +struct generic_product_impl<Lhs, Rhs, SparseShape, SparseTriangularShape, ProductType> + : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType> +{}; + +// sparse-triangular * sparse +template<typename Lhs, typename Rhs, int ProductType> +struct generic_product_impl<Lhs, Rhs, SparseTriangularShape, SparseShape, ProductType> + : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType> +{}; + template<typename Lhs, typename Rhs, int Options> struct evaluator<SparseView<Product<Lhs, Rhs, Options> > > : public evaluator<typename Product<Lhs, Rhs, DefaultProduct>::PlainObject>::type
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index fa9be54..366e272 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp
@@ -194,7 +194,7 @@ VERIFY_IS_APPROX(d3=d1*m2.transpose(), refM3=d1*refM2.transpose()); } - // test self adjoint products + // test self-adjoint and traingular-view products { DenseMatrix b = DenseMatrix::Random(rows, rows); DenseMatrix x = DenseMatrix::Random(rows, rows); @@ -202,9 +202,12 @@ DenseMatrix refUp = DenseMatrix::Zero(rows, rows); DenseMatrix refLo = DenseMatrix::Zero(rows, rows); DenseMatrix refS = DenseMatrix::Zero(rows, rows); + DenseMatrix refA = DenseMatrix::Zero(rows, rows); SparseMatrixType mUp(rows, rows); SparseMatrixType mLo(rows, rows); SparseMatrixType mS(rows, rows); + SparseMatrixType mA(rows, rows); + initSparse<Scalar>(density, refA, mA); do { initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular); } while (refUp.isZero()); @@ -224,19 +227,30 @@ VERIFY_IS_APPROX(mS, refS); VERIFY_IS_APPROX(x=mS*b, refX=refS*b); + // sparse selfadjointView with dense matrices VERIFY_IS_APPROX(x=mUp.template selfadjointView<Upper>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mLo.template selfadjointView<Lower>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mS.template selfadjointView<Upper|Lower>()*b, refX=refS*b); - // sparse selfadjointView * sparse + // sparse selfadjointView with sparse matrices SparseMatrixType mSres(rows,rows); VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>()*mS, refX = refLo.template selfadjointView<Lower>()*refS); - // sparse * sparse selfadjointview VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView<Lower>(), refX = refS * refLo.template selfadjointView<Lower>()); + + // sparse triangularView with dense matrices + VERIFY_IS_APPROX(x=mA.template triangularView<Upper>()*b, refX=refA.template triangularView<Upper>()*b); + VERIFY_IS_APPROX(x=mA.template triangularView<Lower>()*b, refX=refA.template triangularView<Lower>()*b); + VERIFY_IS_APPROX(x=b*mA.template triangularView<Upper>(), refX=b*refA.template triangularView<Upper>()); + VERIFY_IS_APPROX(x=b*mA.template triangularView<Lower>(), refX=b*refA.template triangularView<Lower>()); + + // sparse triangularView with sparse matrices + VERIFY_IS_APPROX(mSres = mA.template triangularView<Lower>()*mS, refX = refA.template triangularView<Lower>()*refS); + VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Lower>(), refX = refS * refA.template triangularView<Lower>()); + VERIFY_IS_APPROX(mSres = mA.template triangularView<Upper>()*mS, refX = refA.template triangularView<Upper>()*refS); + VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Upper>(), refX = refS * refA.template triangularView<Upper>()); } - } // New test for Bug in SparseTimeDenseProduct