Dear all,

I followed the wiki to implement a custom scalar type with support for operator+, operator-, and operator*. Addition and subtraction work fine. However, multiplication yields an error during compilation which I do not know how to fix.

Code: Select all

#include <iostream>
# include <Eigen/Dense>

namespace fdbb {

template<typename T, std::size_t n, std::size_t m>
struct ScalarType
{
    T data[n][m];

    friend std::ostream& operator<< (std::ostream& stream, const ScalarType<T,n,m>& obj) {
        for (std::size_t i=0; i<n; i++)
            for (std::size_t j=0; j<m; j++)
                stream << obj.data[i][j] << " ";
        return stream;
    }
};

template<typename T, std::size_t n, std::size_t m>
inline ScalarType<T,n,m> operator+(const ScalarType<T,n,m>& x,
                                   const ScalarType<T,n,m>& y)
{
    ScalarType<T,n,m> temp;
    for (std::size_t i=0; i<n; i++)
        for (std::size_t j=0; j<m; j++)
            temp.data[i][j] = x.data[i][j] + y.data[i][j];
    return temp;
}

template<typename T, std::size_t n, std::size_t m>
inline ScalarType<T,n,m> operator-(const ScalarType<T,n,m>& x,
                                   const ScalarType<T,n,m>& y)
{
    ScalarType<T,n,m> temp;
    for (std::size_t i=0; i<n; i++)
        for (std::size_t j=0; j<m; j++)
            temp.data[i][j] = x.data[i][j] - y.data[i][j];
    return temp;
}

template<typename T, std::size_t n, std::size_t m, std::size_t o>
inline ScalarType<T,n,o> operator*(const ScalarType<T,n,m>& x,
                                   const ScalarType<T,m,o>& y)
{
    ScalarType<T,n,o> temp;
    for (std::size_t i=0; i<n; i++)
        for (std::size_t j=0; j<o; j++)
            for (std::size_t l=0; l<m; l++)
                temp.data[i][j] += x.data[i][l] * y.data[l][j];
    return temp;
}

} // namespace fdbb

namespace Eigen {

template<typename T, std::size_t n, std::size_t m> struct NumTraits<fdbb::ScalarType<T,n,m>>
 : NumTraits<T> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
  typedef fdbb::ScalarType<T,n,m> Real;
  typedef fdbb::ScalarType<T,n,m> NonInteger;
  typedef fdbb::ScalarType<T,n,m> Nested;
  enum {
    IsComplex = 0,
    IsInteger = 0,
    IsSigned  = 1,
    RequireInitialization = 1,
    ReadCost =   n*m,
    AddCost  = 3*n*m,
    MulCost  = 3*n*m
  };
};

} // namespace Eigen

int main()
{
    constexpr int n = 1;
    constexpr int m = 1;

    typedef fdbb::ScalarType<double,3,3> matrix_type;
   
    // Create 1x1 matrix with 3x3 custom scalar type
    Eigen::Matrix<matrix_type,n,m> A;
    A(0,0) = matrix_type{{{1,2,3},{4,5,6},{7,8,9}}};
   
    std::cout << "Matrix: " << A << std::endl;
    std::cout << "A+A: " << A+A << std::endl;
    std::cout << "A-A: " << A-A << std::endl;
    std::cout << "A*A: " << A*A << std::endl;

    return 0;
}


Compiled under macOS (10.11.6) with

Code: Select all

c++ -g -fpic -std=c++11 -Wall -Wno-unused-function -Wno-attributes -O2 -funroll-loops -ftree-vectorize -ffast-math -pthread -DNDEBUG -D_REENTRANT -D__BLAS -D__EIGEN -D__EIGEN_UNSUPPORTED -I. -I../../include -I/Users/mmoller/codes/libxsmm/samples/eigen/eigen -march=core-avx2 -c /Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp -o build/eigen_test-cpp.o


and

Code: Select all

c++ --version
Apple LLVM version 8.0.0 (clang-800.0.42.1)
Target: x86_64-apple-darwin15.6.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin


I end up with the following error

Code: Select all

In file included from /Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp:10:
In file included from /Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/Dense:1:
In file included from /Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/Core:457:
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/Redux.h:452:12: error: no matching conversion for functional-style cast from 'int' to 'Scalar' (aka 'fdbb::ScalarType<double, 3, 3>')
    return Scalar(0);
           ^~~~~~~~
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/ProductEvaluators.h:251:61: note: in instantiation of member function 'Eigen::DenseBase<Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<fdbb::ScalarType<double, 3,
      3>, fdbb::ScalarType<double, 3, 3> >, const Eigen::Transpose<const Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1> >, const Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1> > >::sum' requested here
    dst.coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
                                                            ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/ProductEvaluators.h:148:37: note: in instantiation of function template specialization 'Eigen::internal::generic_product_impl<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1,
      1, 0, 1, 1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, Eigen::DenseShape, Eigen::DenseShape, 6>::evalTo<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1> >' requested here
    generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
                                    ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/AssignEvaluator.h:836:46: note: in instantiation of member function 'Eigen::internal::Assignment<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>,
      Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, 0>, Eigen::internal::assign_op<fdbb::ScalarType<double, 3, 3>, fdbb::ScalarType<double, 3, 3>
      >, Eigen::internal::Dense2Dense, void>::run' requested here
  Assignment<ActualDstTypeCleaned,Src,Func>::run(actualDst, src, func);
                                             ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/PlainObjectBase.h:728:17: note: in instantiation of function template specialization 'Eigen::internal::call_assignment_no_alias<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>,
      1, 1, 0, 1, 1>, Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, 0>, Eigen::internal::assign_op<fdbb::ScalarType<double, 3, 3>,
      fdbb::ScalarType<double, 3, 3> > >' requested here
      internal::call_assignment_no_alias(this->derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
                ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/PlainObjectBase.h:812:13: note: in instantiation of function template specialization 'Eigen::PlainObjectBase<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>
      >::_set_noalias<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, 0> >' requested here
      this->_set_noalias(other);
            ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/Matrix.h:296:22: note: in instantiation of function template specialization 'Eigen::PlainObjectBase<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>
      >::_init1<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, 0>, Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>,
      Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, 0> >' requested here
      Base::template _init1<T>(x);
                     ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/DenseBase.h:406:14: note: in instantiation of function template specialization 'Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1,
      1>::Matrix<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, 0> >' requested here
      return typename internal::eval<Derived>::type(derived());
             ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/IO.h:220:38: note: in instantiation of member function 'Eigen::DenseBase<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>,
      Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, 0> >::eval' requested here
  return internal::print_matrix(s, m.eval(), EIGEN_DEFAULT_IO_FORMAT);
                                     ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp:103:26: note: in instantiation of function template specialization 'Eigen::operator<<<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>,
      Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, 1, 1, 0, 1, 1>, 0> >' requested here
    std::cout << "A*A: " << A*A << std::endl;
                         ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp:16:8: note: candidate constructor (the implicit copy constructor) not viable: no known conversion from 'int' to 'const fdbb::ScalarType<double, 3, 3>' for 1st argument
struct ScalarType
       ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp:16:8: note: candidate constructor (the implicit move constructor) not viable: no known conversion from 'int' to 'fdbb::ScalarType<double, 3, 3>' for 1st argument
/Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp:20:5: note: candidate constructor not viable: requires 0 arguments, but 1 was provided
    ScalarType() = default;
    ^


Any help is highly appreciated.

Best regards,
Matthias

You need to add a ctor from T, e.g.:

Code: Select all

ScalarType(const T& x) {
   Array<T,n,m>::Map(data) = x;
}


this is needed to handle cases such as:

Matrix<...>(3,0) * Matrix<...>(0,5)

but since in your case sizes are known at compile-time we could remove this code at compile-time, perhaps for the next 3.3.5 release.

Thank you very much for your reply.

I have added the ctor as suggested with

Code: Select all

Eigen::Array<T,n,m>::Map(&data[0][0]) = x;


instead of

Code: Select all

Eigen::Array<T,n,m>::Map(data) = x;


The latter produced the compiler error: no matching member function for call to 'Map'
Matrix-matrix multiplication works well, thanks.

The only downside of this approach is that the default ctors are no longer available by default.
I have implemented

Code: Select all

ScalarType() = default;
ScalarType(const ScalarType<T,n,m>& other) = default;
ScalarType(ScalarType<T,n,m>&& other) = default;
   
ScalarType(std::initializer_list<T> l) {
    std::copy(l.begin(), l.end(), &data[0][0]);
};

ScalarType<T,n,m>& operator=(const ScalarType<T,n,m>& other) = default;


but this only supports matrix construction via

Code: Select all

A(0,0) = matrix_type{1,2,3,4,5,6,7,8,9};


and not the more readable

Code: Select all

A(0,0) = matrix_type{{{1,2,3},{4,5,6},{7,8,9}}};


Any suggestions are appreciated.

Best, Matthias

In addition to the above mentioned problems I encountered the following problem when combining custom scalar types of different type.
Using the implementation of ScalarType from my previous post, the following code in main

Code: Select all

typedef fdbb::ScalarType<double,3,3> matrix_type;
typedef fdbb::ScalarType<double,3,1> vector_type;
   
Eigen::Matrix<matrix_type,1,1> A;
Eigen::Matrix<vector_type,1,1> b;
   
A(0,0) = matrix_type{1,2,3,4,5,6,7,8,9};
b(0)   = vector_type{1,2,3};


leads to the compiler error:

Code: Select all

/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/Product.h:29:116: error: no type named 'ReturnType' in 'Eigen::ScalarBinaryOpTraits<fdbb::ScalarType<double, 3, 3>, fdbb::ScalarType<double, 3, 1>,
      Eigen::internal::scalar_product_op<fdbb::ScalarType<double, 3, 3>, fdbb::ScalarType<double, 3, 1> > >'
  typedef typename ScalarBinaryOpTraits<typename traits<LhsCleaned>::Scalar, typename traits<RhsCleaned>::Scalar>::ReturnType Scalar;


Matrix-vector multiplication of the custom scalar type does work:

Code: Select all

std::cout << matrix_type{1,2,3,4,5,6,7,8,9} * vector_type{1,2,3} << std::endl;


Best,
Matthias

You need to tell Eigen that they are compatible and what the return type is by specializing Eigen::ScalarBinaryOpTraits: http://eigen.tuxfamily.org/dox/structEi ... raits.html

If we move to C++14 and forget about C++03 support, then such type deduction could be done automatically.

Thanks very much. I have implemented the trait as given in the wiki and it seems to work.
My apologise for reporting an error, which turned out to be just a typo in my implementation.

Best, Matthias

I still have some problems with making custom scalar types run. The complete code is given below.

Compilation with the following settings (1x1 dense matrices with compile-time sizes, lines 280-282)

Code: Select all

c++ -g -fpic -std=c++11 -Wall -Wno-unused-function -Wno-attributes -O2 -funroll-loops -ftree-vectorize -ffast-math -pthread -DNDEBUG -D_REENTRANT -D__EIGEN -D__EIGEN_M=1 -D__EIGEN_N=1 -D__EIGEN_K=1 -I. -I../../include -I/Users/mmoller/codes/libxsmm/samples/eigen/eigen -march=core-avx2 -c /Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp -o build/eigen_test-cpp.o


works fine. However, when I switch to dynamic matrix sizes (-D__EIGEN_DYNAMIC, lines 284-286)

Code: Select all

c++ -g -fpic -std=c++11 -Wall -Wno-unused-function -Wno-attributes -O2 -funroll-loops -ftree-vectorize -ffast-math -pthread -DNDEBUG -D_REENTRANT -D__EIGEN -D__EIGEN_M=1 -D__EIGEN_N=1 -D__EIGEN_K=1 -D__EIGEN_DYNAMIC -I. -I../../include -I/Users/mmoller/codes/libxsmm/samples/eigen/eigen -march=core-avx2 -c /Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp -o build/eigen_test-cpp.o


then I get the following compiler errors:

Code: Select all

/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/products/GeneralMatrixMatrix.h:467:32: error: invalid operands to binary expression ('const Scalar'
      (aka 'const ScalarType<double, 3, 1>') and 'const Scalar' (aka 'const fdbb::ScalarType<double, 3, 3>'))
    Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(a_lhs)
                         ~~~~~ ^ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/products/GeneralMatrixMatrix.h:435:7: note: in instantiation of function template specialization
      'Eigen::internal::generic_product_impl<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, -1, -1, 0, -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1,
      -1>, Eigen::DenseShape, Eigen::DenseShape, 8>::scaleAndAddTo<Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1> >' requested here
      scaleAndAddTo(dst, lhs, rhs, Scalar(1));
      ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/ProductEvaluators.h:148:37: note: in instantiation of function template specialization
      'Eigen::internal::generic_product_impl<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, -1, -1, 0, -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1,
      -1>, Eigen::DenseShape, Eigen::DenseShape, 8>::evalTo<Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1> >' requested here
    generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
                                    ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/AssignEvaluator.h:836:46: note: in instantiation of member function
      'Eigen::internal::Assignment<Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, -1, -1, 0,
      -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, 0>, Eigen::internal::assign_op<fdbb::ScalarType<double, 3, 1>, fdbb::ScalarType<double, 3,
      1> >, Eigen::internal::Dense2Dense, void>::run' requested here
  Assignment<ActualDstTypeCleaned,Src,Func>::run(actualDst, src, func);
                                             ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/PlainObjectBase.h:728:17: note: in instantiation of function template specialization
      'Eigen::internal::call_assignment_no_alias<Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3,
      3>, -1, -1, 0, -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, 0>, Eigen::internal::assign_op<fdbb::ScalarType<double, 3, 1>,
      fdbb::ScalarType<double, 3, 1> > >' requested here
      internal::call_assignment_no_alias(this->derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::Scalar>());
                ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/PlainObjectBase.h:812:13: note: in instantiation of function template specialization
      'Eigen::PlainObjectBase<Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1> >::_set_noalias<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>,
      -1, -1, 0, -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, 0> >' requested here
      this->_set_noalias(other);
            ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/Matrix.h:296:22: note: (skipping 1 context in backtrace; use -ftemplate-backtrace-limit=0 to see all)
      Base::template _init1<T>(x);
                     ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/AssignEvaluator.h:796:41: note: in instantiation of function template specialization
      'Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>::Matrix<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, -1, -1, 0, -1, -1>,
      Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, 0> >' requested here
  typename plain_matrix_type<Src>::type tmp(src);
                                        ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/AssignEvaluator.h:782:3: note: in instantiation of function template specialization
      'Eigen::internal::call_assignment<Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, -1,
      -1, 0, -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, 0>, Eigen::internal::assign_op<fdbb::ScalarType<double, 3, 1>,
      fdbb::ScalarType<double, 3, 1> > >' requested here
  call_assignment(dst, src, internal::assign_op<typename Dst::Scalar,typename Src::Scalar>());
  ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/PlainObjectBase.h:710:17: note: in instantiation of function template specialization
      'Eigen::internal::call_assignment<Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, -1,
      -1, 0, -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, 0> >' requested here
      internal::call_assignment(this->derived(), other.derived());
                ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen/Eigen/src/Core/Matrix.h:225:20: note: in instantiation of function template specialization
      'Eigen::PlainObjectBase<Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1> >::_set<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, -1, -1,
      0, -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0, -1, -1>, 0> >' requested here
      return Base::_set(other);
                   ^
/Users/mmoller/codes/libxsmm/samples/eigen/eigen_test.cpp:308:11: note: in instantiation of function template specialization 'Eigen::Matrix<fdbb::ScalarType<double, 3, 1>,
      -1, -1, 0, -1, -1>::operator=<Eigen::Product<Eigen::Matrix<fdbb::ScalarType<double, 3, 3>, -1, -1, 0, -1, -1>, Eigen::Matrix<fdbb::ScalarType<double, 3, 1>, -1, -1, 0,
      -1, -1>, 0> >' requested here
        C = A*B;


Any help is appreciated. Please note that this error is not due to libxsmm which is deactivated in the aforementioned compiler configuration.

Code:

Code: Select all

#include <iostream>
#include <initializer_list>
#include <libxsmm_source.h>

#if !defined(__EIGEN)
/*# define __EIGEN*/
#endif

#if defined(__EIGEN)
# include <Eigen/Dense>
#endif

namespace fdbb {

/**
 * Libxsmm matrix-matrix multipier
 */
template<typename T, int m, int n, int k>
struct ScalarTypeXSMM
{
    static libxsmm_mmfunction<T> xmm;
    static void init();
};

// Initialization of matrix-matrix multiply
template<typename T, int m, int n, int k>
libxsmm_mmfunction<T> ScalarTypeXSMM<T,m,n,k>::xmm;

// Initialization of matrix-matrix multiply
template<typename T, int m, int n, int k>
void ScalarTypeXSMM<T,m,n,k>::init()
{
    xmm = libxsmm_mmfunction<T>(LIBXSMM_GEMM_FLAG_NONE, m, n, k,
                                1/*alpha*/, 0/*beta*/, LIBXSMM_PREFETCH_NONE);
}

/**
 * Custom scalar type
 */
template<typename T, int rows, int cols>
struct ScalarType
{
    // Matrix data
    T data[rows][cols];

    // Default constructor
    ScalarType() = default;

    // Copy constructor
    ScalarType(const ScalarType<T,rows,cols>& other) = default;
   
    // Initializer list constructor
    ScalarType(std::initializer_list<T> l) {
        std::copy(l.begin(), l.end(), &data[0][0]);
    };

    // Scalar constructor (needed by Eigen)
    ScalarType(const T& x) {
        Eigen::Array<T,rows,cols>::Map(&data[0][0]) = x;
    }

    // Copy assignment operator
    ScalarType<T,rows,cols>& operator=(const ScalarType<T,rows,cols>& other) = default;

    // Print to string
    friend std::ostream& operator<< (std::ostream& stream, const ScalarType<T,rows,cols>& obj) {
        stream << "[";
        for (int i=0; i<rows; i++) {
            for (int j=0; j<cols; j++) {
                stream << obj.data[i][j] << (j==cols-1 ? (i==rows-1 ? "" : ";") : ",");
            }
        }
        stream << "]";       
        return stream;
    }

    // Add-assignment operator
    inline ScalarType<T,rows,cols>& operator+=(const ScalarType<T,rows,cols>& x)
    {
        for (int i=0; i<rows; i++)
            for (int j=0; j<cols; j++)
                data[i][j] += x.data[i][j];
        return *this;
    }

    // Sub-assignment operator
    inline ScalarType<T,rows,cols>& operator-=(const ScalarType<T,rows,cols>& x)
    {
        for (int i=0; i<rows; i++)
            for (int j=0; j<cols; j++)
                data[i][j] -= x.data[i][j];
        return *this;
    }
};

// Addition of two custom scalar types
template<typename T, int rows, int cols>
inline ScalarType<T,rows,cols> operator+(const ScalarType<T,rows,cols>& x,
                                         const ScalarType<T,rows,cols>& y)
{
    return ScalarType<T,rows,cols>(x) += y;
}

// Addition of a custom scalar type and a scalar type
template<typename T, int rows, int cols>
inline ScalarType<T,rows,cols> operator+(const ScalarType<T,rows,cols>& x,
                                         const T& y)
{
    ScalarType<T,rows,cols> temp;
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            temp.data[i][j] = x.data[i][j] + y;
    return temp;
}

// Addition of a scalar type and a custom scalar type
template<typename T, int rows, int cols>
inline ScalarType<T,rows,cols> operator+(const T& x,
                                         const ScalarType<T,rows,cols>& y)
{
    ScalarType<T,rows,cols> temp;
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            temp.data[i][j] = x + y.data[i][j];
    return temp;
}

// Subtraction of two custom scalar types
template<typename T, int rows, int cols>
inline ScalarType<T,rows,cols> operator-(const ScalarType<T,rows,cols>& x,
                                         const ScalarType<T,rows,cols>& y)
{
    return ScalarType<T,rows,cols>(x) -= y;
}

// Subtraction of a custom scalar type and a scalar type
template<typename T, int rows, int cols>
inline ScalarType<T,rows,cols> operator-(const ScalarType<T,rows,cols>& x,
                                         const T& y)
{
    ScalarType<T,rows,cols> temp;
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            temp.data[i][j] = x.data[i][j] - y;
    return temp;
}

// Subtraction of a scalar type and a custom scalar type
template<typename T, int rows, int cols>
inline ScalarType<T,rows,cols> operator-(const T& x,
                                         const ScalarType<T,rows,cols>& y)
{
    ScalarType<T,rows,cols> temp;
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            temp.data[i][j] = x - y.data[i][j];
    return temp;
}

// Multiplication of two custom scalar types
template<typename T, int rows, int cols, int k>
inline ScalarType<T,rows,cols> operator*(const ScalarType<T,rows,k>& x,
                                         const ScalarType<T,k,cols>& y)
{
#if defined(__XSMM)
    ScalarType<T,rows,cols> temp;
    ScalarTypeXSMM<T,cols,rows,k>::xmm(&y.data[0][0], &x.data[0][0], &temp.data[0][0]);
#else
    ScalarType<T,rows,cols> temp(0);
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            for (int l=0; l<k; l++)
                temp.data[i][j] += x.data[i][l] * y.data[l][j];
#endif
    return temp;
}

// Multiplication of a custom scalar type and a scalar type
template<typename T, int rows, int cols>
inline ScalarType<T,rows,cols> operator*(const ScalarType<T,rows,cols>& x,
                                         const T& y)
{
    ScalarType<T,rows,cols> temp;
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            temp.data[i][j] = x.data[i][j] * y;
    return temp;
}

// Multiplication of a scalar type and a custom scalar type
template<typename T, int rows, int cols>
inline ScalarType<T,rows,cols> operator*(const T& x,
                                         const ScalarType<T,rows,cols>& y)
{
    ScalarType<T,rows,cols> temp;
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            temp.data[i][j] = x * y.data[i][j];
    return temp;
}

} // namespace fdbb

namespace Eigen {

template<typename T, int rows, int cols>
struct NumTraits<fdbb::ScalarType<T,rows,cols> >
 : NumTraits<T> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
  typedef fdbb::ScalarType<T,rows,cols> Real;
  typedef fdbb::ScalarType<T,rows,cols> NonInteger;
  typedef fdbb::ScalarType<T,rows,cols> Nested;
  enum {
    IsComplex = 0,
    IsInteger = 0,
    IsSigned  = 1,
    RequireInitialization = 1,
    ReadCost =   rows*cols,
    AddCost  = 3*rows*cols,
    MulCost  = 3*rows*cols
  };
};

// Return type specialization for multiplication of two scalar custom types
template<typename T, int rows, int k, int cols>
struct ScalarBinaryOpTraits<typename std::enable_if<rows!=k || k!=cols,fdbb::ScalarType<T,rows,k> >::type,
                            typename std::enable_if<rows!=k || k!=cols,fdbb::ScalarType<T,k,cols> >::type,
                            internal::scalar_product_op<fdbb::ScalarType<T,rows,k>,
                                                        fdbb::ScalarType<T,k,cols> > >
{
    typedef fdbb::ScalarType<T,rows,cols> ReturnType;
};

// Return type specialization for addition of two scalar custom types
template<typename T, int rows, int cols>
struct ScalarBinaryOpTraits<fdbb::ScalarType<T,rows,cols>,
                            fdbb::ScalarType<T,rows,cols>,
                            internal::scalar_sum_op<fdbb::ScalarType<T,rows,cols>,
                                                    fdbb::ScalarType<T,rows,cols> > >
{
    typedef fdbb::ScalarType<T,rows,cols> ReturnType;
};

} // namespace Eigen

int main()
{
    // initialize LIBXSMM
    libxsmm_init();

#define m 3
#define k 3
#define n 1
   
    typedef fdbb::ScalarType<double,m,k> matrix_type;
    typedef fdbb::ScalarType<double,k,n> vector_type;
    typedef fdbb::ScalarType<double,m,n> result_type;
   
    fdbb::ScalarTypeXSMM<double,n,m,k>::init();

#ifndef __EIGEN_M
#define M 1
#else
#define M __EIGEN_M
#endif

#ifndef __EIGEN_N
#define N 1
#else
#define N __EIGEN_N
#endif

#ifndef __EIGEN_K
#define K 1
#else
#define K __EIGEN_K
#endif

#if !defined(__EIGEN_DYNAMIC)
    Eigen::Matrix<matrix_type,M,K> A;
    Eigen::Matrix<vector_type,K,N> B;
    Eigen::Matrix<result_type,M,N> C;
#else
    Eigen::Matrix<matrix_type,Eigen::Dynamic,Eigen::Dynamic> A(M,K);
    Eigen::Matrix<vector_type,Eigen::Dynamic,Eigen::Dynamic> B(K,N);
    Eigen::Matrix<result_type,Eigen::Dynamic,Eigen::Dynamic> C(M,N);
#endif
   
    for (int i=0; i<M; i++) {
        for (int j=0; j<K; j++) {
            A(i,j) = double(i+j+1) * matrix_type{1,2,3,4,5,6,7,8,9};
        }
    }
   
    for (int i=0; i<K; i++) {
        for (int j=0; j<N; j++) {
            B(i,j) = double(i+j+1) * vector_type{1,2,3,4,5,6,7,8,9};
        }
    }
   
    std::cout << "A = \n" << A << std::endl;

    std::cout << "B = \n" << B << std::endl;
   
    const unsigned long long start = libxsmm_timer_tick();

    for (int run=0; run<8000000; run++)
        C = A*B;
       
    const unsigned long long end = libxsmm_timer_tick();
    const double duration = libxsmm_timer_duration(start, end);

    std::cout << "Duration: " << 1000.0 * duration << std::endl;

    std::cout << "C = \n" << C << std::endl;
   
    // finalize LIBXSMM
    libxsmm_finalize();
   
    return 0;
}


use A.lazyProduct(B) to by pass the difficulties.

Thanks very much. Since lazy evaluation did not work as expected I had a deeper look into the Eigen kernel

IMHO, there is a bug in Eigen/src/SparseCore/ConservativeSparseSparseProduct.h.

For custom scalar types it is not correct to use the same Scalar type for the Lhs, Rhs, and the ResultType in the functions conservative_sparse_sparse_product_impl and sparse_sparse_to_dense_product_impl. I have attached a patch file that makes the necessary changes:

PATCH FILE

Code: Select all

diff -ruN eigen-eigen-5a0156e40feb/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h eigen/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h
--- eigen-eigen-5a0156e40feb/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h   2017-06-15 09:10:20.000000000 +0200
+++ eigen/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h   2017-08-31 12:39:50.000000000 +0200
@@ -17,7 +17,9 @@
 template<typename Lhs, typename Rhs, typename ResultType>
 static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res, bool sortedInsertion = false)
 {
-  typedef typename remove_all<Lhs>::type::Scalar Scalar;
+  typedef typename remove_all<Lhs>::type::Scalar LhsScalar;
+  typedef typename remove_all<Rhs>::type::Scalar RhsScalar;
+  typedef typename remove_all<ResultType>::type::Scalar Scalar;
 
   // make sure to call innerSize/outerSize since we fake the storage order.
   Index rows = lhs.innerSize();
@@ -51,12 +53,12 @@
     Index nnz = 0;
     for (typename evaluator<Rhs>::InnerIterator rhsIt(rhsEval, j); rhsIt; ++rhsIt)
     {
-      Scalar y = rhsIt.value();
+      RhsScalar y = rhsIt.value();
       Index k = rhsIt.index();
       for (typename evaluator<Lhs>::InnerIterator lhsIt(lhsEval, k); lhsIt; ++lhsIt)
       {
         Index i = lhsIt.index();
-        Scalar x = lhsIt.value();
+        LhsScalar x = lhsIt.value();
         if(!mask[i])
         {
           mask[i] = true;
@@ -263,7 +265,8 @@
 template<typename Lhs, typename Rhs, typename ResultType>
 static void sparse_sparse_to_dense_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
 {
-  typedef typename remove_all<Lhs>::type::Scalar Scalar;
+  typedef typename remove_all<Lhs>::type::Scalar LhsScalar;
+  typedef typename remove_all<Rhs>::type::Scalar RhsScalar;
   Index cols = rhs.outerSize();
   eigen_assert(lhs.outerSize() == rhs.innerSize());
 
@@ -274,12 +277,12 @@
   {
     for (typename evaluator<Rhs>::InnerIterator rhsIt(rhsEval, j); rhsIt; ++rhsIt)
     {
-      Scalar y = rhsIt.value();
+      RhsScalar y = rhsIt.value();
       Index k = rhsIt.index();
       for (typename evaluator<Lhs>::InnerIterator lhsIt(lhsEval, k); lhsIt; ++lhsIt)
       {
         Index i = lhsIt.index();
-        Scalar x = lhsIt.value();
+        LhsScalar x = lhsIt.value();
         res.coeffRef(i,j) += x * y;
       }


thank you for the feedback, there were more occurrences though: https://bitbucket.org/eigen/eigen/commits/bccf2d933699