Re-enable products with triangular views of sparse matrices: we simply have to treat them as a sparse matrix.

This commit is contained in:
Gael Guennebaud 2014-10-06 16:11:26 +02:00
parent 893bfcf95f
commit d44d432baa
5 changed files with 49 additions and 8 deletions

View File

@ -247,7 +247,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
inline const AdjointReturnType adjoint() const inline const AdjointReturnType adjoint() const
{ return AdjointReturnType(m_matrix.adjoint()); } { return AdjointReturnType(m_matrix.adjoint()); }
typedef TriangularView<Transpose<MatrixType>,TransposeMode> TransposeReturnType; typedef TriangularView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType;
/** \sa MatrixBase::transpose() */ /** \sa MatrixBase::transpose() */
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
inline TransposeReturnType transpose() inline TransposeReturnType transpose()
@ -255,11 +255,13 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
EIGEN_STATIC_ASSERT_LVALUE(MatrixType) EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
return TransposeReturnType(m_matrix.const_cast_derived().transpose()); return TransposeReturnType(m_matrix.const_cast_derived().transpose());
} }
typedef TriangularView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType;
/** \sa MatrixBase::transpose() const */ /** \sa MatrixBase::transpose() const */
EIGEN_DEVICE_FUNC 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> template<typename Other>

View File

@ -141,6 +141,8 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,C
typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrixAux; typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrixAux;
typedef typename sparse_eval<ColMajorMatrixAux,ResultType::RowsAtCompileTime,ResultType::ColsAtCompileTime>::type ColMajorMatrix; 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. // FIXME, the following heuristic is probably not very good.
if(lhs.rows()>=rhs.cols()) if(lhs.rows()>=rhs.cols())
{ {

View File

@ -146,6 +146,11 @@ struct generic_product_impl<Lhs, Rhs, SparseShape, DenseShape, ProductType>
} }
}; };
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> template<typename Lhs, typename Rhs, int ProductType>
struct generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType> struct generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
{ {
@ -158,12 +163,17 @@ struct generic_product_impl<Lhs, Rhs, DenseShape, SparseShape, ProductType>
RhsNested rhsNested(rhs); RhsNested rhsNested(rhs);
dst.setZero(); dst.setZero();
// transpoe everything // transpose everything
Transpose<Dest> dstT(dst); Transpose<Dest> dstT(dst);
internal::sparse_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, typename Dest::Scalar(1)); 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> template<typename LhsT, typename RhsT, bool NeedToTranspose>
struct sparse_dense_outer_product_evaluator struct sparse_dense_outer_product_evaluator
{ {

View File

@ -33,6 +33,7 @@ SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other
namespace internal { namespace internal {
// sparse * sparse
template<typename Lhs, typename Rhs, int ProductType> template<typename Lhs, typename Rhs, int ProductType>
struct generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType> struct generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
{ {
@ -48,6 +49,18 @@ struct generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
} }
}; };
// 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> template<typename Lhs, typename Rhs, int Options>
struct evaluator<SparseView<Product<Lhs, Rhs, Options> > > struct evaluator<SparseView<Product<Lhs, Rhs, Options> > >
: public evaluator<typename Product<Lhs, Rhs, DefaultProduct>::PlainObject>::type : public evaluator<typename Product<Lhs, Rhs, DefaultProduct>::PlainObject>::type

View File

@ -194,7 +194,7 @@ template<typename SparseMatrixType> void sparse_product()
VERIFY_IS_APPROX(d3=d1*m2.transpose(), refM3=d1*refM2.transpose()); 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 b = DenseMatrix::Random(rows, rows);
DenseMatrix x = DenseMatrix::Random(rows, rows); DenseMatrix x = DenseMatrix::Random(rows, rows);
@ -202,9 +202,12 @@ template<typename SparseMatrixType> void sparse_product()
DenseMatrix refUp = DenseMatrix::Zero(rows, rows); DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
DenseMatrix refLo = DenseMatrix::Zero(rows, rows); DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
DenseMatrix refS = DenseMatrix::Zero(rows, rows); DenseMatrix refS = DenseMatrix::Zero(rows, rows);
DenseMatrix refA = DenseMatrix::Zero(rows, rows);
SparseMatrixType mUp(rows, rows); SparseMatrixType mUp(rows, rows);
SparseMatrixType mLo(rows, rows); SparseMatrixType mLo(rows, rows);
SparseMatrixType mS(rows, rows); SparseMatrixType mS(rows, rows);
SparseMatrixType mA(rows, rows);
initSparse<Scalar>(density, refA, mA);
do { do {
initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular); initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
} while (refUp.isZero()); } while (refUp.isZero());
@ -224,19 +227,30 @@ template<typename SparseMatrixType> void sparse_product()
VERIFY_IS_APPROX(mS, refS); VERIFY_IS_APPROX(mS, refS);
VERIFY_IS_APPROX(x=mS*b, refX=refS*b); 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=mUp.template selfadjointView<Upper>()*b, refX=refS*b);
VERIFY_IS_APPROX(x=mLo.template selfadjointView<Lower>()*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); 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); SparseMatrixType mSres(rows,rows);
VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>()*mS, VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>()*mS,
refX = refLo.template selfadjointView<Lower>()*refS); refX = refLo.template selfadjointView<Lower>()*refS);
// sparse * sparse selfadjointview
VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView<Lower>(), VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView<Lower>(),
refX = refS * refLo.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 // New test for Bug in SparseTimeDenseProduct