ufunc: use Eigen algebra to compute dotfunc

Author: jcarpent

include/eigenpy/user-type.hpp

@@ -199,18 +199,15 @@ namespace eigenpy
                                  void * op, npy_intp n, void * /*arr*/)
       {
 //        std::cout << "dotfunc" << std::endl;
-        T res(0);
-        char *ip0 = (char*)ip0_, *ip1 = (char*)ip1_;
-        npy_intp i;
-        for(i = 0; i < n; i++)
-        {
-          
-          res += *static_cast<T*>(static_cast<void*>(ip0))
-          * *static_cast<T*>(static_cast<void*>(ip1));
-          ip0 += is0;
-          ip1 += is1;
-        }
-        *static_cast<T*>(op) = res;
+        typedef Eigen::Matrix<T,Eigen::Dynamic,1> VectorT;
+        typedef Eigen::InnerStride<Eigen::Dynamic> InputStride;
+        typedef const Eigen::Map<const VectorT,0,InputStride> ConstMapType;
+
+        ConstMapType
+        v0(static_cast<T*>(ip0_),n,InputStride(is0/sizeof(T))),
+        v1(static_cast<T*>(ip1_),n,InputStride(is1/sizeof(T)));
+        
+        *static_cast<T*>(op) = v0.dot(v1);
       }
 
 

unittest/python/test_user_type.py

@@ -7,6 +7,7 @@
 
 def test(dtype):
   mat = np.ones((rows,cols),dtype=dtype)
+  mat = np.random.rand(rows,cols).astype(dtype)
   mat_copy = mat.copy()
   assert (mat == mat_copy).all()
   assert not (mat != mat_copy).all()
@@ -33,8 +34,11 @@ def test(dtype):
   assert not (mat < mat).all()
 
   mat2 = mat.dot(mat.T)
+  mat2_ref = mat.astype(np.double).dot(mat.T.astype(np.double))
+  assert np.isclose(mat2.astype(np.double),mat2_ref).all()
   if np.__version__ >= '1.17.0':
     mat2 = np.matmul(mat,mat.T)
+    assert np.isclose(mat2.astype(np.double),mat2_ref).all()
 
 def test_cast(from_dtype,to_dtype):
   np.can_cast(from_dtype,to_dtype)

unittest/user_type.cpp

@@ -200,6 +200,8 @@ BOOST_PYTHON_MODULE(user_type)
   eigenpy::registerCast<int32_t,DoubleType>(true);
   eigenpy::registerCast<DoubleType,int64_t>(false);
   eigenpy::registerCast<int64_t,DoubleType>(true);
+  eigenpy::registerCast<FloatType,double>(true);
+  eigenpy::registerCast<double,FloatType>(false);
   eigenpy::registerCast<FloatType,int64_t>(false);
   eigenpy::registerCast<int64_t,FloatType>(true);