Fix Map-with-Stride and cover it by new unit tests.

Eigen/src/Core/DenseStorageBase.h

@@ -362,6 +362,9 @@ class DenseStorageBase : public _Base<Derived>
       * while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned
       * \a data pointers.
       *
+      * These methods do not allow to specify strides. If you need to specify strides, you have to
+      * use the Map class directly.
+      *
       * \see class Map
       */
     //@{

Eigen/src/Core/Map.h

@@ -52,11 +52,23 @@ template<typename MatrixType, int Options, typename StrideType>
 struct ei_traits<Map<MatrixType, Options, StrideType> >
   : public ei_traits<MatrixType>
 {
+  typedef typename MatrixType::Scalar Scalar;
   enum {
+    InnerStride = StrideType::InnerStrideAtCompileTime,
+    OuterStride = StrideType::OuterStrideAtCompileTime,
+    HasNoInnerStride = InnerStride <= 1,
+    HasNoOuterStride = OuterStride == 0,
+    HasNoStride = HasNoInnerStride && HasNoOuterStride,
+    IsAligned = int(int(Options)&Aligned)==Aligned,
+    IsDynamicSize = MatrixType::SizeAtCompileTime==Dynamic,
+    KeepsPacketAccess = bool(HasNoInnerStride)
+                        && ( bool(IsDynamicSize)
+                           || HasNoOuterStride
+                           || ( OuterStride!=Dynamic && ((int(OuterStride)*sizeof(Scalar))%16)==0 ) ),
     Flags0 = ei_traits<MatrixType>::Flags,
-    Flags1 = ((Options&Aligned)==Aligned ? Flags0 |  AlignedBit
-                                         : Flags0 & ~AlignedBit),
-    Flags = int(StrideType::InnerStrideAtCompileTime)==1 ? Flags1 : (Flags1 & ~PacketAccessBit)
+    Flags1 = IsAligned ? int(Flags0) |  AlignedBit : int(Flags0) & ~AlignedBit,
+    Flags2 = HasNoStride ? int(Flags1) : int(Flags1 & ~LinearAccessBit),
+    Flags = KeepsPacketAccess ? int(Flags2) : (int(Flags2) & ~PacketAccessBit)
   };
 };
 
@@ -94,23 +106,6 @@ template<typename MatrixType, int Options, typename StrideType> class Map
     inline Map(const Scalar* data, int rows, int cols, const StrideType& stride = StrideType())
       : Base(data, rows, cols), m_stride(stride) {}
 
-/*
-    inline void resize(int rows, int cols)
-    {
-      EIGEN_ONLY_USED_FOR_DEBUG(rows);
-      EIGEN_ONLY_USED_FOR_DEBUG(cols);
-      ei_assert(rows == this->rows());
-      ei_assert(cols == this->cols());
-    }
-
-    inline void resize(int size)
-    {
-      EIGEN_STATIC_ASSERT_VECTOR_ONLY(MatrixType)
-      EIGEN_ONLY_USED_FOR_DEBUG(size);
-      ei_assert(size == this->size());
-    }
-*/
-
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Map)
 
   protected:

Eigen/src/Core/Stride.h

@@ -25,7 +25,7 @@
 #ifndef EIGEN_STRIDE_H
 #define EIGEN_STRIDE_H
 
-template<int _InnerStrideAtCompileTime, int _OuterStrideAtCompileTime>
+template<int _OuterStrideAtCompileTime, int _InnerStrideAtCompileTime>
 class Stride
 {
   public:
@@ -36,45 +36,45 @@ class Stride
     };
 
     Stride()
-      : m_inner(InnerStrideAtCompileTime), m_outer(OuterStrideAtCompileTime)
+      : m_outer(OuterStrideAtCompileTime), m_inner(InnerStrideAtCompileTime)
     {
       ei_assert(InnerStrideAtCompileTime != Dynamic && OuterStrideAtCompileTime != Dynamic);
     }
 
-    Stride(int innerStride, int outerStride)
-      : m_inner(innerStride), m_outer(outerStride)
+    Stride(int outerStride, int innerStride)
+      : m_outer(outerStride), m_inner(innerStride)
     {
       ei_assert(innerStride>=0 && outerStride>=0);
     }
 
     Stride(const Stride& other)
-      : m_inner(other.inner()), m_outer(other.outer())
+      : m_outer(other.outer()), m_inner(other.inner())
     {}
 
-    inline int inner() const { return m_inner.value(); }
     inline int outer() const { return m_outer.value(); }
+    inline int inner() const { return m_inner.value(); }
 
   protected:
-    ei_int_if_dynamic<InnerStrideAtCompileTime> m_inner;
     ei_int_if_dynamic<OuterStrideAtCompileTime> m_outer;
+    ei_int_if_dynamic<InnerStrideAtCompileTime> m_inner;
 };
 
-template<int Value = Dynamic>
-class InnerStride : public Stride<Value, 0>
-{
-    typedef Stride<Value,0> Base;
-  public:
-    InnerStride() : Base() {}
-    InnerStride(int v) : Base(v,0) {}
-};
-
-template<int Value = Dynamic>
-class OuterStride : public Stride<0, Value>
+template<int Value>
+class InnerStride : public Stride<0, Value>
 {
     typedef Stride<0, Value> Base;
+  public:
+    InnerStride() : Base() {}
+    InnerStride(int v) : Base(0, v) {}
+};
+
+template<int Value>
+class OuterStride : public Stride<Value, 0>
+{
+    typedef Stride<Value, 0> Base;
   public:
     OuterStride() : Base() {}
-    OuterStride(int v) : Base(0,v) {}
+    OuterStride(int v) : Base(v,0) {}
 };
 
 #endif // EIGEN_STRIDE_H

Eigen/src/Core/util/Constants.h

@@ -86,11 +86,11 @@ const unsigned int EvalBeforeAssigningBit = 0x4;
   * Long version: means that the coefficients can be handled by packets
   * and start at a memory location whose alignment meets the requirements
   * of the present CPU architecture for optimized packet access. In the fixed-size
-  * case, there is the additional condition that the total size of the coefficients
-  * array is a multiple of the packet size, so that it is possible to access all the
-  * coefficients by packets. In the dynamic-size case, there is no such condition
-  * on the total size, so it might not be possible to access the few last coeffs
-  * by packets.
+  * case, there is the additional condition that it be possible to access all the
+  * coefficients by packets (this implies the requirement that the size be a multiple of 16 bytes,
+  * and that any nontrivial strides don't break the alignment). In the dynamic-size case,
+  * there is no such condition on the total size and strides, so it might not be possible to access
+  * all coeffs by packets.
   *
   * \note This bit can be set regardless of whether vectorization is actually enabled.
   *       To check for actual vectorizability, see \a ActualPacketAccessBit.

Eigen/src/Core/util/ForwardDeclarations.h

@@ -61,7 +61,7 @@ template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime=SizeA
 template<typename MatrixType, typename DiagonalType, int ProductOrder> class DiagonalProduct;
 template<typename MatrixType, int Index> class Diagonal;
 
-template<int InnerStrideAtCompileTime = Dynamic, int OuterStrideAtCompileTime = Dynamic> class Stride;
+template<int InnerStrideAtCompileTime, int OuterStrideAtCompileTime> class Stride;
 template<typename MatrixType, int Options=Unaligned, typename StrideType = Stride<0,0> > class Map;
 
 template<typename Derived> class TriangularBase;

test/CMakeLists.txt

@@ -115,6 +115,7 @@ ei_add_test(miscmatrices)
 ei_add_test(commainitializer)
 ei_add_test(smallvectors)
 ei_add_test(map)
+ei_add_test(mapstride)
 ei_add_test(array)
 ei_add_test(array_for_matrix)
 ei_add_test(array_replicate)

test/map.cpp

@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2006-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -42,8 +42,8 @@ template<typename VectorType> void map_class_vector(const VectorType& m)
   VectorType ma1 = Map<VectorType, Aligned>(array1, size);
   VectorType ma2 = Map<VectorType, Aligned>(array2, size);
   VectorType ma3 = Map<VectorType>(array3unaligned, size);
-  VERIFY_IS_APPROX(ma1, ma2);
-  VERIFY_IS_APPROX(ma1, ma3);
+  VERIFY_IS_EQUAL(ma1, ma2);
+  VERIFY_IS_EQUAL(ma1, ma3);
   VERIFY_RAISES_ASSERT((Map<VectorType,Aligned>(array3unaligned, size)));
 
   ei_aligned_delete(array1, size);
@@ -70,9 +70,9 @@ template<typename MatrixType> void map_class_matrix(const MatrixType& m)
   Map<MatrixType>(array3unaligned, rows, cols) = Map<MatrixType>(array1, rows, cols);
   MatrixType ma1 = Map<MatrixType>(array1, rows, cols);
   MatrixType ma2 = Map<MatrixType, Aligned>(array2, rows, cols);
-  VERIFY_IS_APPROX(ma1, ma2);
+  VERIFY_IS_EQUAL(ma1, ma2);
   MatrixType ma3 = Map<MatrixType>(array3unaligned, rows, cols);
-  VERIFY_IS_APPROX(ma1, ma3);
+  VERIFY_IS_EQUAL(ma1, ma3);
 
   ei_aligned_delete(array1, size);
   ei_aligned_delete(array2, size);
@@ -97,8 +97,8 @@ template<typename VectorType> void map_static_methods(const VectorType& m)
   VectorType ma1 = VectorType::Map(array1, size);
   VectorType ma2 = VectorType::MapAligned(array2, size);
   VectorType ma3 = VectorType::Map(array3unaligned, size);
-  VERIFY_IS_APPROX(ma1, ma2);
-  VERIFY_IS_APPROX(ma1, ma3);
+  VERIFY_IS_EQUAL(ma1, ma2);
+  VERIFY_IS_EQUAL(ma1, ma3);
 
   ei_aligned_delete(array1, size);
   ei_aligned_delete(array2, size);

test/mapstride.cpp

@@ -0,0 +1,139 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+template<typename VectorType> void map_class_vector(const VectorType& m)
+{
+  typedef typename VectorType::Scalar Scalar;
+
+  int size = m.size();
+
+  VectorType v = VectorType::Random(size);
+
+  int arraysize = 3*size;
+  
+  Scalar* array = ei_aligned_new<Scalar>(arraysize);
+
+  {
+    Map<VectorType, Aligned, InnerStride<3> > map(array, size);
+    map = v;
+    for(int i = 0; i < size; ++i)
+    {
+      VERIFY(array[3*i] == v[i]);
+      VERIFY(map[i] == v[i]);
+    }
+  }
+
+  {
+    Map<VectorType, Unaligned, InnerStride<Dynamic> > map(array, size, InnerStride<Dynamic>(2));
+    map = v;
+    for(int i = 0; i < size; ++i)
+    {
+      VERIFY(array[2*i] == v[i]);
+      VERIFY(map[i] == v[i]);
+    }
+  }
+
+  ei_aligned_delete(array, arraysize);
+}
+
+template<typename MatrixType> void map_class_matrix(const MatrixType& _m)
+{
+  typedef typename MatrixType::Scalar Scalar;
+
+  int rows = _m.rows(), cols = _m.cols();
+
+  MatrixType m = MatrixType::Random(rows,cols);
+
+  int arraysize = 2*(rows+4)*(cols+4);
+
+  Scalar* array = ei_aligned_new<Scalar>(arraysize);
+
+  // test no inner stride and some dynamic outer stride
+  {
+    Map<MatrixType, Aligned, OuterStride<Dynamic> > map(array, rows, cols, OuterStride<Dynamic>(m.innerSize()+1));
+    map = m;
+    VERIFY(map.outerStride() == map.innerSize()+1);
+    for(int i = 0; i < m.outerSize(); ++i)
+      for(int j = 0; j < m.innerSize(); ++j)
+      {
+        VERIFY(array[map.outerStride()*i+j] == m.coeffByOuterInner(i,j));
+        VERIFY(map.coeffByOuterInner(i,j) == m.coeffByOuterInner(i,j));
+      }
+  }
+
+  // test no inner stride and an outer stride of +4. This is quite important as for fixed-size matrices,
+  // this allows to hit the special case where it's vectorizable.
+  {
+    enum {
+      InnerSize = MatrixType::InnerSizeAtCompileTime,
+      OuterStrideAtCompileTime = InnerSize==Dynamic ? Dynamic : InnerSize+4
+    };
+    Map<MatrixType, Aligned, OuterStride<OuterStrideAtCompileTime> >
+      map(array, rows, cols, OuterStride<OuterStrideAtCompileTime>(m.innerSize()+4));
+    map = m;
+    VERIFY(map.outerStride() == map.innerSize()+4);
+    for(int i = 0; i < m.outerSize(); ++i)
+      for(int j = 0; j < m.innerSize(); ++j)
+      {
+        VERIFY(array[map.outerStride()*i+j] == m.coeffByOuterInner(i,j));
+        VERIFY(map.coeffByOuterInner(i,j) == m.coeffByOuterInner(i,j));
+      }
+  }
+
+  // test both inner stride and outer stride
+  {
+    Map<MatrixType, Aligned, Stride<Dynamic,Dynamic> > map(array, rows, cols, Stride<Dynamic,Dynamic>(2*m.innerSize()+1, 2));
+    map = m;
+    VERIFY(map.outerStride() == 2*map.innerSize()+1);
+    VERIFY(map.innerStride() == 2);
+    for(int i = 0; i < m.outerSize(); ++i)
+      for(int j = 0; j < m.innerSize(); ++j)
+      {
+        VERIFY(array[map.outerStride()*i+map.innerStride()*j] == m.coeffByOuterInner(i,j));
+        VERIFY(map.coeffByOuterInner(i,j) == m.coeffByOuterInner(i,j));
+      }
+  }
+
+  ei_aligned_delete(array, arraysize);
+}
+
+void test_mapstride()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_1( map_class_vector(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST_2( map_class_vector(Vector4d()) );
+    CALL_SUBTEST_3( map_class_vector(RowVector4f()) );
+    CALL_SUBTEST_4( map_class_vector(VectorXcf(8)) );
+    CALL_SUBTEST_5( map_class_vector(VectorXi(12)) );
+
+    CALL_SUBTEST_1( map_class_matrix(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST_2( map_class_matrix(Matrix4d()) );
+    CALL_SUBTEST_3( map_class_matrix(Matrix<float,3,5>()) );
+    CALL_SUBTEST_3( map_class_matrix(Matrix<float,4,8>()) );
+    CALL_SUBTEST_4( map_class_matrix(MatrixXcf(ei_random<int>(1,10),ei_random<int>(1,10))) );
+    CALL_SUBTEST_5( map_class_matrix(MatrixXi(5,5)));//ei_random<int>(1,10),ei_random<int>(1,10))) );
+  }
+}

test/vectorization_logic.cpp

@@ -33,6 +33,14 @@ bool test_assign(const Dst&, const Src&, int traversal, int unrolling)
     && ei_assign_traits<Dst,Src>::Unrolling==unrolling;
 }
 
+template<typename Dst, typename Src>
+bool test_assign(int traversal, int unrolling)
+{
+  ei_assign_traits<Dst,Src>::debug();
+  return ei_assign_traits<Dst,Src>::Traversal==traversal
+    && ei_assign_traits<Dst,Src>::Unrolling==unrolling;
+}
+
 template<typename Xpr>
 bool test_redux(const Xpr&, int traversal, int unrolling)
 {
@@ -86,6 +94,15 @@ void test_vectorization_logic()
   VERIFY(test_assign(MatrixXf(10,10),MatrixXf(20,20).block(10,10,2,3),
     SliceVectorizedTraversal,NoUnrolling));
 
+  VERIFY((test_assign<
+           Map<Matrix<float,4,8>, Aligned, OuterStride<12> >,
+           Matrix<float,4,8>
+          >(InnerVectorizedTraversal,CompleteUnrolling)));
+
+  VERIFY((test_assign<
+           Map<Matrix<float,4,8>, Aligned, InnerStride<12> >,
+           Matrix<float,4,8>
+          >(DefaultTraversal,CompleteUnrolling)));
 
   VERIFY(test_redux(VectorXf(10),
     LinearVectorizedTraversal,NoUnrolling));