From a16599751f42242a3cbb80a00cddc983a6bb2675 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Mon, 31 Aug 2009 17:39:56 +0200
Subject: fix Matrix::stride for vectors, add a unit test for Block::stride and
 make use of it where it was relevant

---
 Eigen/src/Core/Matrix.h                           | 17 +++++++--
 Eigen/src/Core/products/SelfadjointMatrixVector.h | 12 +++---
 test/submatrices.cpp                              | 45 +++++++++++++++++++++++
 3 files changed, 64 insertions(+), 10 deletions(-)

diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index f58424ba2..d0603871a 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -142,12 +142,21 @@ class Matrix
     EIGEN_STRONG_INLINE int rows() const { return m_storage.rows(); }
     EIGEN_STRONG_INLINE int cols() const { return m_storage.cols(); }
 
-    EIGEN_STRONG_INLINE int stride(void) const
+    /** Returns the leading dimension (for matrices) or the increment (for vectors) to be used with data().
+      *
+      * More precisely:
+      *  - for a column major matrix it returns the number of elements between two successive columns
+      *  - for a row major matrix it returns the number of elements between two successive rows
+      *  - for a vector it returns the number of elements between two successive coefficients
+      * This function has to be used together with the MapBase::data() function.
+      *
+      * \sa Matrix::data() */
+    EIGEN_STRONG_INLINE int stride() const
     {
-      if(Flags & RowMajorBit)
-        return m_storage.cols();
+      if(IsVectorAtCompileTime)
+        return 1;
       else
-        return m_storage.rows();
+        return (Flags & RowMajorBit) ? m_storage.cols() : m_storage.rows();
     }
 
     EIGEN_STRONG_INLINE const Scalar& coeff(int row, int col) const
diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h
index c2c33d5b8..d5927307d 100644
--- a/Eigen/src/Core/products/SelfadjointMatrixVector.h
+++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h
@@ -185,14 +185,14 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
     Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
                                * RhsBlasTraits::extractScalarFactor(m_rhs);
 
-    ei_assert((&dst.coeff(1))-(&dst.coeff(0))==1 && "not implemented yet");
+    ei_assert(dst.stride()==1 && "not implemented yet");
     ei_product_selfadjoint_vector<Scalar, ei_traits<_ActualLhsType>::Flags&RowMajorBit, int(LhsUpLo), bool(LhsBlasTraits::NeedToConjugate), bool(RhsBlasTraits::NeedToConjugate)>
       (
-        lhs.rows(),                                     // size
-        &lhs.coeff(0,0), lhs.stride(),                  // lhs info
-        &rhs.coeff(0), (&rhs.coeff(1))-(&rhs.coeff(0)), // rhs info
-        &dst.coeffRef(0),                               // result info
-        actualAlpha                                     // scale factor
+        lhs.rows(),                       // size
+        &lhs.coeff(0,0),  lhs.stride(),   // lhs info
+        &rhs.coeff(0),    rhs.stride(),   // rhs info
+        &dst.coeffRef(0),                 // result info
+        actualAlpha                       // scale factor
       );
   }
 };
diff --git a/test/submatrices.cpp b/test/submatrices.cpp
index a819cadc2..6fe86c281 100644
--- a/test/submatrices.cpp
+++ b/test/submatrices.cpp
@@ -170,6 +170,48 @@ template<typename MatrixType> void submatrices(const MatrixType& m)
   VERIFY(ei_real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols));
 }
 
+
+template<typename MatrixType>
+void compare_using_data_and_stride(const MatrixType& m)
+{
+  int rows = m.rows();
+  int cols = m.cols();
+  int size = m.size();
+  int stride = m.stride();
+  const typename MatrixType::Scalar* data = m.data();
+
+  for(int j=0;j<cols;++j)
+    for(int i=0;i<rows;++i)
+      VERIFY_IS_APPROX(m.coeff(i,j), data[(MatrixType::Flags&RowMajorBit) ? i*stride+j : j*stride + i]);
+
+  if(MatrixType::IsVectorAtCompileTime)
+  {
+    VERIFY_IS_APPROX(stride, int((&m.coeff(1))-(&m.coeff(0))));
+    for (int i=0;i<size;++i)
+      VERIFY_IS_APPROX(m.coeff(i), data[i*stride]);
+  }
+}
+
+template<typename MatrixType>
+void data_and_stride(const MatrixType& m)
+{
+  int rows = m.rows();
+  int cols = m.cols();
+
+  int r1 = ei_random<int>(0,rows-1);
+  int r2 = ei_random<int>(r1,rows-1);
+  int c1 = ei_random<int>(0,cols-1);
+  int c2 = ei_random<int>(c1,cols-1);
+
+  MatrixType m1 = MatrixType::Random(rows, cols);
+  compare_using_data_and_stride(m1.block(r1, c1, r2-r1+1, c2-c1+1));
+  compare_using_data_and_stride(m1.transpose().block(c1, r1, c2-c1+1, r2-r1+1));
+  compare_using_data_and_stride(m1.row(r1));
+  compare_using_data_and_stride(m1.col(c1));
+  compare_using_data_and_stride(m1.row(r1).transpose());
+  compare_using_data_and_stride(m1.col(c1).transpose());
+}
+
 void test_submatrices()
 {
   for(int i = 0; i < g_repeat; i++) {
@@ -179,5 +221,8 @@ void test_submatrices()
     CALL_SUBTEST( submatrices(MatrixXi(8, 12)) );
     CALL_SUBTEST( submatrices(MatrixXcd(20, 20)) );
     CALL_SUBTEST( submatrices(MatrixXf(20, 20)) );
+
+    CALL_SUBTEST( data_and_stride(MatrixXf(ei_random(5,50), ei_random(5,50))) );
+    CALL_SUBTEST( data_and_stride(Matrix<int,Dynamic,Dynamic,RowMajor>(ei_random(5,50), ei_random(5,50))) );
   }
 }
-- 
cgit v1.2.3