Dear all,

While experimenting with custom scalar types I figured out that dense-sparse matrix-matrix multiply dispatches to sparse-dense matrix-matrix multiply with both matrices transposed (and the resulting matrix transposed, too). However, the scalar types of all matrices involved remain unchanged. I have implemented a custom scalar type that represents local dense n x m matrices including all admissible operations like (n x m) = (n x k) * (k x m). If the global container matrices are transposed then this has to be done for the custom scalar types as well. Otherwise, the global multiply is no longer well defined. Is there an internal trait that I can specialize for my custom scalar type that tells Eigen how to treat the scalar type in case the transpose() method is called on the global matrix?

Best regard,
Matthias

yes, currently we assume fields with associative addition and multiplication, which is not your case. This is also a problem for, say, {Sparse}Matrix<Quaternion>. Such a trick is used is several other matrix product cases. The easiest might be to wrap scalar products like "a*b" through a template helper that would swap the factor if needed.

I looked into the file Eigen/src/SparseCore/SparseDenseProduct.h

1) IMHO, it might be better to replace the occurrence of 'const typename Res::Scalar& alpha' by 'const AlphaType& alpha; in the various run(...) methods since (a) AlphaType is a template parameter of the struct anyway, and (b) the method that should be called by the user

Code: Select all

template<typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType>
inline void sparse_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
 {
   sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, AlphaType>::run(lhs, rhs, res, alpha);
}


follows this convention anyway. See attached patch:

2) Would you mind adding a dense_time_sparse_product_impl or, alternatively, add an extra template parameter to sparse_time_dense_product_impl, which let's the user decide whether internally res = lhs * rhs or res = rhs * lhs should be computed? If I understand you correctly, Eigen assumes associativity of addition and multiplication. In my case the latter is not true but I could very easily re-use the sparse_time_dense_product_impl by simply swapping the order lhs * rhs -> rhs * lhs. If this could be controlled by an extra templated bool flag that would be great.

Patch for issue 1)

Code: Select all

From 4c30bda61429a6a0d417b28d3432421e3ebf631b Mon Sep 17 00:00:00 2001
From: Matthias Moeller <m.moller@tudelft.nl>
Date: Mon, 18 Sep 2017 15:32:39 +0200
Subject: [PATCH] Changed type of scalar alpha coefficient in
 sparse_time_dense_product_impl to template parameter

---
 Eigen/src/SparseCore/SparseDenseProduct.h | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h
index 0547db596..21836612a 100644
--- a/Eigen/src/SparseCore/SparseDenseProduct.h
+++ b/Eigen/src/SparseCore/SparseDenseProduct.h
@@ -23,15 +23,15 @@ template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,
          bool ColPerCol = ((DenseRhsType::Flags&RowMajorBit)==0) || DenseRhsType::ColsAtCompileTime==1>
 struct sparse_time_dense_product_impl;
 
-template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
-struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, RowMajor, true>
+template<typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType>
+struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType,AlphaType, RowMajor, true>
 {
   typedef typename internal::remove_all<SparseLhsType>::type Lhs;
   typedef typename internal::remove_all<DenseRhsType>::type Rhs;
   typedef typename internal::remove_all<DenseResType>::type Res;
   typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
   typedef evaluator<Lhs> LhsEval;
-  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
+  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
   {
     LhsEval lhsEval(lhs);
     
@@ -61,7 +61,7 @@ struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, t
     }
   }
   
-  static void processRow(const LhsEval& lhsEval, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha, Index i, Index col)
+  static void processRow(const LhsEval& lhsEval, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha, Index i, Index col)
   {
     typename Res::Scalar tmp(0);
     for(LhsInnerIterator it(lhsEval,i); it ;++it)
@@ -83,7 +83,7 @@ struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, t
 // };
 
 template<typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType>
-struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, AlphaType, ColMajor, true>
+struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType,AlphaType, ColMajor, true>
 {
   typedef typename internal::remove_all<SparseLhsType>::type Lhs;
   typedef typename internal::remove_all<DenseRhsType>::type Rhs;
@@ -105,14 +105,14 @@ struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, A
   }
 };
 
-template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
-struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, RowMajor, false>
+template<typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType>
+struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType,AlphaType, RowMajor, false>
 {
   typedef typename internal::remove_all<SparseLhsType>::type Lhs;
   typedef typename internal::remove_all<DenseRhsType>::type Rhs;
   typedef typename internal::remove_all<DenseResType>::type Res;
   typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
-  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
+  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
   {
     evaluator<Lhs> lhsEval(lhs);
     for(Index j=0; j<lhs.outerSize(); ++j)
@@ -124,14 +124,14 @@ struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, t
   }
 };
 
-template<typename SparseLhsType, typename DenseRhsType, typename DenseResType>
-struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, typename DenseResType::Scalar, ColMajor, false>
+template<typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType>
+struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType,AlphaType, ColMajor, false>
 {
   typedef typename internal::remove_all<SparseLhsType>::type Lhs;
   typedef typename internal::remove_all<DenseRhsType>::type Rhs;
   typedef typename internal::remove_all<DenseResType>::type Res;
   typedef typename evaluator<Lhs>::InnerIterator LhsInnerIterator;
-  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha)
+  static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
   {
     evaluator<Lhs> lhsEval(lhs);
     for(Index j=0; j<lhs.outerSize(); ++j)
@@ -143,7 +143,7 @@ struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, t
   }
 };
 
-template<typename SparseLhsType, typename DenseRhsType, typename DenseResType,typename AlphaType>
+template<typename SparseLhsType, typename DenseRhsType, typename DenseResType, typename AlphaType>
 inline void sparse_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha)
 {
   sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, AlphaType>::run(lhs, rhs, res, alpha);
--
2.14.1