24 files changed, 430 insertions, 227 deletions

diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h
index 9440cebf1..3133aa03a 100644
--- a/Eigen/src/Core/Assign.h
+++ b/Eigen/src/Core/Assign.h
@@ -37,34 +37,37 @@ struct ei_assign_traits
 public:
   enum {
     DstIsAligned = Derived::Flags & AlignedBit,
+    DstHasDirectAccess = Derived::Flags & DirectAccessBit,
     SrcIsAligned = OtherDerived::Flags & AlignedBit,
-    SrcAlignment = DstIsAligned && SrcIsAligned ? Aligned : Unaligned
+    JointAlignment = DstIsAligned && SrcIsAligned ? Aligned : Unaligned
   };
 
 private:
   enum {
-    InnerSize = int(Derived::Flags)&RowMajorBit
-              ? Derived::ColsAtCompileTime
-              : Derived::RowsAtCompileTime,
-    InnerMaxSize = int(Derived::Flags)&RowMajorBit
-              ? Derived::MaxColsAtCompileTime
-              : Derived::MaxRowsAtCompileTime,
-    MaxSizeAtCompileTime = ei_size_at_compile_time<Derived::MaxColsAtCompileTime,Derived::MaxRowsAtCompileTime>::ret,
+    InnerSize = int(Derived::IsVectorAtCompileTime) ? int(Derived::SizeAtCompileTime)
+              : int(Derived::Flags)&RowMajorBit ? int(Derived::ColsAtCompileTime)
+              : int(Derived::RowsAtCompileTime),
+    InnerMaxSize = int(Derived::IsVectorAtCompileTime) ? int(Derived::MaxSizeAtCompileTime)
+              : int(Derived::Flags)&RowMajorBit ? int(Derived::MaxColsAtCompileTime)
+              : int(Derived::MaxRowsAtCompileTime),
+    MaxSizeAtCompileTime = Derived::SizeAtCompileTime,
     PacketSize = ei_packet_traits<typename Derived::Scalar>::size
   };
 
   enum {
-    StorageOrdersAgree = (int(Derived::Flags)&RowMajorBit)==(int(OtherDerived::Flags)&RowMajorBit),
+    LhsIsEffectivelyRowMajor = (Derived::RowsAtCompileTime==1) || (int(Derived::Flags)&RowMajorBit),
+    RhsIsEffectivelyRowMajor = (OtherDerived::RowsAtCompileTime==1) || (int(OtherDerived::Flags)&RowMajorBit),
+    StorageOrdersAgree = (LhsIsEffectivelyRowMajor == RhsIsEffectivelyRowMajor),
     MightVectorize = StorageOrdersAgree
                   && (int(Derived::Flags) & int(OtherDerived::Flags) & ActualPacketAccessBit),
     MayInnerVectorize  = MightVectorize && int(InnerSize)!=Dynamic && int(InnerSize)%int(PacketSize)==0
                        && int(DstIsAligned) && int(SrcIsAligned),
     MayLinearize = StorageOrdersAgree && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit),
-    MayLinearVectorize = MightVectorize && MayLinearize
-                                        && (DstIsAligned || MaxSizeAtCompileTime == Dynamic),
+    MayLinearVectorize = MightVectorize && MayLinearize && DstHasDirectAccess
+                       && (DstIsAligned || MaxSizeAtCompileTime == Dynamic),
       /* If the destination isn't aligned, we have to do runtime checks and we don't unroll,
          so it's only good for large enough sizes. */
-    MaySliceVectorize  = MightVectorize && int(InnerMaxSize)>=3*PacketSize
+    MaySliceVectorize  = MightVectorize && DstHasDirectAccess && int(InnerMaxSize)>=3*PacketSize
       /* slice vectorization can be slow, so we only want it if the slices are big, which is
          indicated by InnerMaxSize rather than InnerSize, think of the case of a dynamic block
          in a fixed-size matrix */
@@ -108,12 +111,13 @@ public:
   {
     EIGEN_DEBUG_VAR(DstIsAligned)
     EIGEN_DEBUG_VAR(SrcIsAligned)
-    EIGEN_DEBUG_VAR(SrcAlignment)
+    EIGEN_DEBUG_VAR(JointAlignment)
     EIGEN_DEBUG_VAR(InnerSize)
     EIGEN_DEBUG_VAR(InnerMaxSize)
     EIGEN_DEBUG_VAR(PacketSize)
     EIGEN_DEBUG_VAR(StorageOrdersAgree)
     EIGEN_DEBUG_VAR(MightVectorize)
+    EIGEN_DEBUG_VAR(MayLinearize)
     EIGEN_DEBUG_VAR(MayInnerVectorize)
     EIGEN_DEBUG_VAR(MayLinearVectorize)
     EIGEN_DEBUG_VAR(MaySliceVectorize)
@@ -137,17 +141,13 @@ template<typename Derived1, typename Derived2, int Index, int Stop>
 struct ei_assign_DefaultTraversal_CompleteUnrolling
 {
   enum {
-    row = int(Derived1::Flags)&RowMajorBit
-        ? Index / int(Derived1::ColsAtCompileTime)
-        : Index % Derived1::RowsAtCompileTime,
-    col = int(Derived1::Flags)&RowMajorBit
-        ? Index % int(Derived1::ColsAtCompileTime)
-        : Index / Derived1::RowsAtCompileTime
+    outer = Index / Derived1::InnerSizeAtCompileTime,
+    inner = Index % Derived1::InnerSizeAtCompileTime
   };
 
   EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src)
   {
-    dst.copyCoeff(row, col, src);
+    dst.copyCoeffByOuterInner(outer, inner, src);
     ei_assign_DefaultTraversal_CompleteUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src);
   }
 };
@@ -161,13 +161,10 @@ struct ei_assign_DefaultTraversal_CompleteUnrolling<Derived1, Derived2, Stop, St
 template<typename Derived1, typename Derived2, int Index, int Stop>
 struct ei_assign_DefaultTraversal_InnerUnrolling
 {
-  EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src, int row_or_col)
+  EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src, int outer)
   {
-    const bool rowMajor = int(Derived1::Flags)&RowMajorBit;
-    const int row = rowMajor ? row_or_col : Index;
-    const int col = rowMajor ? Index : row_or_col;
-    dst.copyCoeff(row, col, src);
-    ei_assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src, row_or_col);
+    dst.copyCoeffByOuterInner(outer, Index, src);
+    ei_assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src, outer);
   }
 };
 
@@ -205,18 +202,14 @@ template<typename Derived1, typename Derived2, int Index, int Stop>
 struct ei_assign_innervec_CompleteUnrolling
 {
   enum {
-    row = int(Derived1::Flags)&RowMajorBit
-        ? Index / int(Derived1::ColsAtCompileTime)
-        : Index % Derived1::RowsAtCompileTime,
-    col = int(Derived1::Flags)&RowMajorBit
-        ? Index % int(Derived1::ColsAtCompileTime)
-        : Index / Derived1::RowsAtCompileTime,
-    SrcAlignment = ei_assign_traits<Derived1,Derived2>::SrcAlignment
+    outer = Index / Derived1::InnerSizeAtCompileTime,
+    inner = Index % Derived1::InnerSizeAtCompileTime,
+    JointAlignment = ei_assign_traits<Derived1,Derived2>::JointAlignment
   };
 
   EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src)
   {
-    dst.template copyPacket<Derived2, Aligned, SrcAlignment>(row, col, src);
+    dst.template copyPacketByOuterInner<Derived2, Aligned, JointAlignment>(outer, inner, src);
     ei_assign_innervec_CompleteUnrolling<Derived1, Derived2,
       Index+ei_packet_traits<typename Derived1::Scalar>::size, Stop>::run(dst, src);
   }
@@ -231,13 +224,11 @@ struct ei_assign_innervec_CompleteUnrolling<Derived1, Derived2, Stop, Stop>
 template<typename Derived1, typename Derived2, int Index, int Stop>
 struct ei_assign_innervec_InnerUnrolling
 {
-  EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src, int row_or_col)
+  EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src, int outer)
   {
-    const int row = int(Derived1::Flags)&RowMajorBit ? row_or_col : Index;
-    const int col = int(Derived1::Flags)&RowMajorBit ? Index : row_or_col;
-    dst.template copyPacket<Derived2, Aligned, Aligned>(row, col, src);
+    dst.template copyPacketByOuterInner<Derived2, Aligned, Aligned>(outer, Index, src);
     ei_assign_innervec_InnerUnrolling<Derived1, Derived2,
-      Index+ei_packet_traits<typename Derived1::Scalar>::size, Stop>::run(dst, src, row_or_col);
+      Index+ei_packet_traits<typename Derived1::Scalar>::size, Stop>::run(dst, src, outer);
   }
 };
 
@@ -267,14 +258,9 @@ struct ei_assign_impl<Derived1, Derived2, DefaultTraversal, NoUnrolling>
   {
     const int innerSize = dst.innerSize();
     const int outerSize = dst.outerSize();
-    for(int j = 0; j < outerSize; ++j)
-      for(int i = 0; i < innerSize; ++i)
-      {
-        if(int(Derived1::Flags)&RowMajorBit)
-          dst.copyCoeff(j, i, src);
-        else
-          dst.copyCoeff(i, j, src);
-      }
+    for(int outer = 0; outer < outerSize; ++outer)
+      for(int inner = 0; inner < innerSize; ++inner)
+        dst.copyCoeffByOuterInner(outer, inner, src);
   }
 };
 
@@ -293,12 +279,10 @@ struct ei_assign_impl<Derived1, Derived2, DefaultTraversal, InnerUnrolling>
 {
   EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src)
   {
-    const bool rowMajor = int(Derived1::Flags)&RowMajorBit;
-    const int innerSize = rowMajor ? Derived1::ColsAtCompileTime : Derived1::RowsAtCompileTime;
     const int outerSize = dst.outerSize();
-    for(int j = 0; j < outerSize; ++j)
-      ei_assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, 0, innerSize>
-        ::run(dst, src, j);
+    for(int outer = 0; outer < outerSize; ++outer)
+      ei_assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, 0, Derived1::InnerSizeAtCompileTime>
+        ::run(dst, src, outer);
   }
 };
 
@@ -339,14 +323,9 @@ struct ei_assign_impl<Derived1, Derived2, InnerVectorizedTraversal, NoUnrolling>
     const int innerSize = dst.innerSize();
     const int outerSize = dst.outerSize();
     const int packetSize = ei_packet_traits<typename Derived1::Scalar>::size;
-    for(int j = 0; j < outerSize; ++j)
-      for(int i = 0; i < innerSize; i+=packetSize)
-      {
-        if(int(Derived1::Flags)&RowMajorBit)
-          dst.template copyPacket<Derived2, Aligned, Aligned>(j, i, src);
-        else
-          dst.template copyPacket<Derived2, Aligned, Aligned>(i, j, src);
-      }
+    for(int outer = 0; outer < outerSize; ++outer)
+      for(int inner = 0; inner < innerSize; inner+=packetSize)
+        dst.template copyPacketByOuterInner<Derived2, Aligned, Aligned>(outer, inner, src);
   }
 };
 
@@ -365,12 +344,10 @@ struct ei_assign_impl<Derived1, Derived2, InnerVectorizedTraversal, InnerUnrolli
 {
   EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src)
   {
-    const bool rowMajor = int(Derived1::Flags)&RowMajorBit;
-    const int innerSize = rowMajor ? Derived1::ColsAtCompileTime : Derived1::RowsAtCompileTime;
     const int outerSize = dst.outerSize();
-    for(int j = 0; j < outerSize; ++j)
-      ei_assign_innervec_InnerUnrolling<Derived1, Derived2, 0, innerSize>
-        ::run(dst, src, j);
+    for(int outer = 0; outer < outerSize; ++outer)
+      ei_assign_innervec_InnerUnrolling<Derived1, Derived2, 0, Derived1::InnerSizeAtCompileTime>
+        ::run(dst, src, outer);
   }
 };
 
@@ -418,7 +395,7 @@ struct ei_assign_impl<Derived1, Derived2, LinearVectorizedTraversal, NoUnrolling
 
     for(int index = alignedStart; index < alignedEnd; index += packetSize)
     {
-      dst.template copyPacket<Derived2, Aligned, ei_assign_traits<Derived1,Derived2>::SrcAlignment>(index, src);
+      dst.template copyPacket<Derived2, Aligned, ei_assign_traits<Derived1,Derived2>::JointAlignment>(index, src);
     }
 
     ei_unaligned_assign_impl<>::run(src,dst,alignedEnd,size);
@@ -452,40 +429,24 @@ struct ei_assign_impl<Derived1, Derived2, SliceVectorizedTraversal, NoUnrolling>
     const int packetAlignedMask = packetSize - 1;
     const int innerSize = dst.innerSize();
     const int outerSize = dst.outerSize();
-    const int alignedStep = (packetSize - dst.stride() % packetSize) & packetAlignedMask;
+    const int alignedStep = (packetSize - dst.outerStride() % packetSize) & packetAlignedMask;
     int alignedStart = ei_assign_traits<Derived1,Derived2>::DstIsAligned ? 0
                      : ei_first_aligned(&dst.coeffRef(0,0), innerSize);
 
-    for(int i = 0; i < outerSize; ++i)
+    for(int outer = 0; outer < outerSize; ++outer)
     {
       const int alignedEnd = alignedStart + ((innerSize-alignedStart) & ~packetAlignedMask);
-
       // do the non-vectorizable part of the assignment
-      for (int index = 0; index<alignedStart ; ++index)
-      {
-        if(Derived1::Flags&RowMajorBit)
-          dst.copyCoeff(i, index, src);
-        else
-          dst.copyCoeff(index, i, src);
-      }
+      for(int inner = 0; inner<alignedStart ; ++inner)
+        dst.copyCoeffByOuterInner(outer, inner, src);
 
       // do the vectorizable part of the assignment
-      for (int index = alignedStart; index<alignedEnd; index+=packetSize)
-      {
-        if(Derived1::Flags&RowMajorBit)
-          dst.template copyPacket<Derived2, Aligned, Unaligned>(i, index, src);
-        else
-          dst.template copyPacket<Derived2, Aligned, Unaligned>(index, i, src);
-      }
+      for(int inner = alignedStart; inner<alignedEnd; inner+=packetSize)
+        dst.template copyPacketByOuterInner<Derived2, Aligned, Unaligned>(outer, inner, src);
 
       // do the non-vectorizable part of the assignment
-      for (int index = alignedEnd; index<innerSize ; ++index)
-      {
-        if(Derived1::Flags&RowMajorBit)
-          dst.copyCoeff(i, index, src);
-        else
-          dst.copyCoeff(index, i, src);
-      }
+      for(int inner = alignedEnd; inner<innerSize ; ++inner)
+        dst.copyCoeffByOuterInner(outer, inner, src);
 
       alignedStart = std::min<int>((alignedStart+alignedStep)%packetSize, innerSize);
     }
diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h
index 3b4234c22..d3c4dfa99 100644
--- a/Eigen/src/Core/Block.h
+++ b/Eigen/src/Core/Block.h
@@ -100,8 +100,8 @@ template<typename MatrixType, int BlockRows, int BlockCols, int _DirectAccessSta
         // The case a 1x1 matrix seems ambiguous, but the result is the same anyway.
         m_startRow( (BlockRows==1) && (BlockCols==MatrixType::ColsAtCompileTime) ? i : 0),
         m_startCol( (BlockRows==MatrixType::RowsAtCompileTime) && (BlockCols==1) ? i : 0),
-        m_blockRows(matrix.rows()), // if it is a row, then m_blockRows has a fixed-size of 1, so no pb to try to overwrite it
-        m_blockCols(matrix.cols())  // same for m_blockCols
+        m_blockRows(BlockRows==1 ? 1 : matrix.rows()),
+        m_blockCols(BlockCols==1 ? 1 : matrix.cols())
     {
       ei_assert( (i>=0) && (
           ((BlockRows==1) && (BlockCols==MatrixType::ColsAtCompileTime) && i<matrix.rows())
@@ -112,7 +112,7 @@ template<typename MatrixType, int BlockRows, int BlockCols, int _DirectAccessSta
       */
     inline Block(const MatrixType& matrix, int startRow, int startCol)
       : m_matrix(matrix), m_startRow(startRow), m_startCol(startCol),
-        m_blockRows(matrix.rows()), m_blockCols(matrix.cols())
+        m_blockRows(BlockRows), m_blockCols(BlockCols)
     {
       EIGEN_STATIC_ASSERT(RowsAtCompileTime!=Dynamic && ColsAtCompileTime!=Dynamic,THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE)
       ei_assert(startRow >= 0 && BlockRows >= 1 && startRow + BlockRows <= matrix.rows()
@@ -196,8 +196,8 @@ template<typename MatrixType, int BlockRows, int BlockCols, int _DirectAccessSta
     #ifdef EIGEN_PARSED_BY_DOXYGEN
     /** \sa MapBase::data() */
     inline const Scalar* data() const;
-    /** \sa MapBase::stride() */
-    inline int stride() const;
+    inline int innerStride() const;
+    inline int outerStride() const;
     #endif
 
   protected:
@@ -260,17 +260,23 @@ class Block<MatrixType,BlockRows,BlockCols,HasDirectAccess>
              && startCol >= 0 && blockCols >= 0 && startCol + blockCols <= matrix.cols());
     }
 
-    /** \sa MapBase::stride() */
-    inline int stride() const
+    /** \sa MapBase::innerStride() */
+    inline int innerStride() const
     {
-      return    ((!Base::IsVectorAtCompileTime)
-              || (BlockRows==1 && ((Flags&RowMajorBit)==0))
-              || (BlockCols==1 && ((Flags&RowMajorBit)==RowMajorBit)))
-              ? m_matrix.stride() : 1;
+      return RowsAtCompileTime==1 ? m_matrix.colStride()
+           : ColsAtCompileTime==1 ? m_matrix.rowStride()
+           : m_matrix.innerStride();
+    }
+    
+    /** \sa MapBase::outerStride() */
+    inline int outerStride() const
+    {
+      return IsVectorAtCompileTime ? this->size() : m_matrix.outerStride();
     }
 
   #ifndef __SUNPRO_CC
   // FIXME sunstudio is not friendly with the above friend...
+  // META-FIXME there is no 'friend' keyword around here. Is this obsolete?
   protected:
   #endif
 
diff --git a/Eigen/src/Core/Coeffs.h b/Eigen/src/Core/Coeffs.h
index ebfd0c80e..da7b9153f 100644
--- a/Eigen/src/Core/Coeffs.h
+++ b/Eigen/src/Core/Coeffs.h
@@ -25,6 +25,24 @@
 #ifndef EIGEN_COEFFS_H
 #define EIGEN_COEFFS_H
 
+template<typename Derived>
+EIGEN_STRONG_INLINE int DenseBase<Derived>::rowIndexByOuterInner(int outer, int inner)
+{
+  return int(Derived::RowsAtCompileTime) == 1 ? 0
+       : int(Derived::ColsAtCompileTime) == 1 ? inner
+       : int(Derived::Flags)&RowMajorBit ? outer
+       : inner;
+}
+
+template<typename Derived>
+EIGEN_STRONG_INLINE int DenseBase<Derived>::colIndexByOuterInner(int outer, int inner)
+{
+  return int(Derived::ColsAtCompileTime) == 1 ? 0
+       : int(Derived::RowsAtCompileTime) == 1 ? inner
+       : int(Derived::Flags)&RowMajorBit ? inner
+       : outer;
+}
+
 /** Short version: don't use this function, use
   * \link operator()(int,int) const \endlink instead.
   *
@@ -48,6 +66,14 @@ EIGEN_STRONG_INLINE const typename DenseBase<Derived>::CoeffReturnType DenseBase
   return derived().coeff(row, col);
 }
 
+template<typename Derived>
+EIGEN_STRONG_INLINE const typename DenseBase<Derived>::CoeffReturnType DenseBase<Derived>
+  ::coeffByOuterInner(int outer, int inner) const
+{
+  return coeff(rowIndexByOuterInner(outer, inner),
+               colIndexByOuterInner(outer, inner));
+}
+
 /** \returns the coefficient at given the given row and column.
   *
   * \sa operator()(int,int), operator[](int) const
@@ -84,6 +110,14 @@ EIGEN_STRONG_INLINE typename ei_traits<Derived>::Scalar& DenseBase<Derived>
   return derived().coeffRef(row, col);
 }
 
+template<typename Derived>
+EIGEN_STRONG_INLINE typename ei_traits<Derived>::Scalar& DenseBase<Derived>
+  ::coeffRefByOuterInner(int outer, int inner)
+{
+  return coeffRef(rowIndexByOuterInner(outer, inner),
+                  colIndexByOuterInner(outer, inner));
+}
+
 /** \returns a reference to the coefficient at given the given row and column.
   *
   * \sa operator()(int,int) const, operator[](int)
@@ -261,6 +295,15 @@ DenseBase<Derived>::packet(int row, int col) const
   return derived().template packet<LoadMode>(row,col);
 }
 
+template<typename Derived>
+template<int LoadMode>
+EIGEN_STRONG_INLINE typename ei_packet_traits<typename ei_traits<Derived>::Scalar>::type
+DenseBase<Derived>::packetByOuterInner(int outer, int inner) const
+{
+  return packet<LoadMode>(rowIndexByOuterInner(outer, inner),
+                          colIndexByOuterInner(outer, inner));
+}
+
 /** Stores the given packet of coefficients, at the given row and column of this expression. It is your responsibility
   * to ensure that a packet really starts there. This method is only available on expressions having the
   * PacketAccessBit.
@@ -279,6 +322,16 @@ EIGEN_STRONG_INLINE void DenseBase<Derived>::writePacket
   derived().template writePacket<StoreMode>(row,col,x);
 }
 
+template<typename Derived>
+template<int StoreMode>
+EIGEN_STRONG_INLINE void DenseBase<Derived>::writePacketByOuterInner
+(int outer, int inner, const typename ei_packet_traits<typename ei_traits<Derived>::Scalar>::type& x)
+{
+  writePacket<StoreMode>(rowIndexByOuterInner(outer, inner),
+                         colIndexByOuterInner(outer, inner),
+                         x);
+}
+
 /** \returns the packet of coefficients starting at the given index. It is your responsibility
   * to ensure that a packet really starts there. This method is only available on expressions having the
   * PacketAccessBit and the LinearAccessBit.
@@ -346,6 +399,16 @@ EIGEN_STRONG_INLINE void DenseBase<Derived>::copyCoeff(int index, const DenseBas
   derived().coeffRef(index) = other.derived().coeff(index);
 }
 
+template<typename Derived>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE void DenseBase<Derived>::copyCoeffByOuterInner(int outer, int inner, const DenseBase<OtherDerived>& other)
+{
+  const int row = Derived::rowIndexByOuterInner(outer,inner);
+  const int col = Derived::colIndexByOuterInner(outer,inner);
+  // derived() is important here: copyCoeff() may be reimplemented in Derived!
+  derived().copyCoeff(row, col, other);
+}
+
 /** \internal Copies the packet at position (row,col) of other into *this.
   *
   * This method is overridden in SwapWrapper, allowing swap() assignments to share 99% of their code
@@ -379,6 +442,16 @@ EIGEN_STRONG_INLINE void DenseBase<Derived>::copyPacket(int index, const DenseBa
     other.derived().template packet<LoadMode>(index));
 }
 
+template<typename Derived>
+template<typename OtherDerived, int StoreMode, int LoadMode>
+EIGEN_STRONG_INLINE void DenseBase<Derived>::copyPacketByOuterInner(int outer, int inner, const DenseBase<OtherDerived>& other)
+{
+  const int row = Derived::rowIndexByOuterInner(outer,inner);
+  const int col = Derived::colIndexByOuterInner(outer,inner);
+  // derived() is important here: copyCoeff() may be reimplemented in Derived!
+  derived().copyPacket<OtherDerived, StoreMode, LoadMode>(row, col, other);
+}
+
 template<typename Derived, bool JustReturnZero>
 struct ei_first_aligned_impl
 {
diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h
index eb018738f..5682d7278 100644
--- a/Eigen/src/Core/DenseBase.h
+++ b/Eigen/src/Core/DenseBase.h
@@ -126,6 +126,9 @@ template<typename Derived> class DenseBase
 
       IsRowMajor = int(Flags) & RowMajorBit, /**< True if this expression is row major. */
 
+      InnerSizeAtCompileTime = int(IsVectorAtCompileTime) ? SizeAtCompileTime
+                             : int(Flags)&RowMajorBit ? ColsAtCompileTime : RowsAtCompileTime,
+
       CoeffReadCost = ei_traits<Derived>::CoeffReadCost,
         /**< This is a rough measure of how expensive it is to read one coefficient from
           * this expression.
@@ -160,31 +163,104 @@ template<typename Derived> class DenseBase
       * In other words, this function returns
       * \code rows()==1 || cols()==1 \endcode
       * \sa rows(), cols(), IsVectorAtCompileTime. */
-    inline bool isVector() const { return rows()==1 || cols()==1; }
-    /** \returns the size of the storage major dimension,
-      * i.e., the number of columns for a columns major matrix, and the number of rows otherwise */
-    int outerSize() const { return (int(Flags)&RowMajorBit) ? this->rows() : this->cols(); }
-    /** \returns the size of the inner dimension according to the storage order,
-      * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */
-    int innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); }
-
-    /** Only plain matrices, not expressions may be resized; therefore the only useful resize method is
-      * Matrix::resize(). The present method only asserts that the new size equals the old size, and does
+
+    /** \returns the outer size.
+      *
+      * \note For a vector, this returns just 1. For a matrix (non-vector), this is the major dimension
+      * with respect to the storage order, i.e., the number of columns for a column-major matrix,
+      * and the number of rows for a row-major matrix. */
+    int outerSize() const
+    {
+      return IsVectorAtCompileTime ? 1
+           : (int(Flags)&RowMajorBit) ? this->rows() : this->cols();
+    }
+
+    /** \returns the inner size.
+      *
+      * \note For a vector, this is just the size. For a matrix (non-vector), this is the minor dimension
+      * with respect to the storage order, i.e., the number of rows for a column-major matrix,
+      * and the number of columns for a row-major matrix. */
+    int innerSize() const
+    {
+      return IsVectorAtCompileTime ? this->size()
+           : (int(Flags)&RowMajorBit) ? this->cols() : this->rows();
+    }
+
+    /** Only plain matrices/arrays, not expressions, may be resized; therefore the only useful resize methods are
+      * Matrix::resize() and Array::resize(). The present method only asserts that the new size equals the old size, and does
       * nothing else.
       */
     void resize(int size)
     {
       ei_assert(size == this->size()
-                && "MatrixBase::resize() does not actually allow to resize.");
+                && "DenseBase::resize() does not actually allow to resize.");
     }
-    /** Only plain matrices, not expressions may be resized; therefore the only useful resize method is
-      * Matrix::resize(). The present method only asserts that the new size equals the old size, and does
+    /** Only plain matrices/arrays, not expressions, may be resized; therefore the only useful resize methods are
+      * Matrix::resize() and Array::resize(). The present method only asserts that the new size equals the old size, and does
       * nothing else.
       */
     void resize(int rows, int cols)
     {
       ei_assert(rows == this->rows() && cols == this->cols()
-                && "MatrixBase::resize() does not actually allow to resize.");
+                && "DenseBase::resize() does not actually allow to resize.");
+    }
+
+    /** \returns the pointer increment between two consecutive elements.
+      *
+      * \note For vectors, the storage order is ignored. For matrices (non-vectors), we're looking
+      *       at the increment between two consecutive elements within a slice in the inner direction.
+      *
+      * \sa outerStride(), rowStride(), colStride()
+      */
+    inline int innerStride() const
+    {
+      EIGEN_STATIC_ASSERT(int(Flags)&DirectAccessBit,
+                          THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
+      return derived().innerStride();
+    }
+
+    /** \returns the pointer increment between two consecutive inner slices (for example, between two consecutive columns
+      *          in a column-major matrix).
+      *
+      * \note For vectors, the storage order is ignored, there is only one inner slice, and so this method returns 1.
+      *       For matrices (non-vectors), the notion of inner slice depends on the storage order.
+      *
+      * \sa innerStride(), rowStride(), colStride()
+      */
+    inline int outerStride() const
+    {
+      EIGEN_STATIC_ASSERT(int(Flags)&DirectAccessBit,
+                          THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
+      return derived().outerStride();
+    }
+
+    inline int stride() const
+    {
+      return IsVectorAtCompileTime ? innerStride() : outerStride();
+    }
+
+    /** \returns the pointer increment between two consecutive rows.
+      *
+      * \sa innerStride(), outerStride(), colStride()
+      */
+    inline int rowStride() const
+    {
+      return ColsAtCompileTime==1 ? innerStride()
+           : RowsAtCompileTime==1 ? outerStride()
+           : IsRowMajor ? outerStride()
+           : innerStride();
+    }
+
+    /** \returns the pointer increment between two consecutive columns.
+      *
+      * \sa innerStride(), outerStride(), rowStride()
+      */
+    inline int colStride() const
+    {
+      return ColsAtCompileTime==1 ? outerStride()
+           : RowsAtCompileTime==1 ? innerStride()
+           : IsRowMajor ? innerStride()
+           : outerStride();
     }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
@@ -242,9 +318,11 @@ template<typename Derived> class DenseBase
     CommaInitializer<Derived> operator<< (const DenseBase<OtherDerived>& other);
 
     const CoeffReturnType coeff(int row, int col) const;
+    const CoeffReturnType coeffByOuterInner(int outer, int inner) const;
     const CoeffReturnType operator()(int row, int col) const;
 
     Scalar& coeffRef(int row, int col);
+    Scalar& coeffRefByOuterInner(int outer, int inner);
     Scalar& operator()(int row, int col);
 
     const CoeffReturnType coeff(int index) const;
@@ -259,17 +337,30 @@ template<typename Derived> class DenseBase
     template<typename OtherDerived>
     void copyCoeff(int row, int col, const DenseBase<OtherDerived>& other);
     template<typename OtherDerived>
+    void copyCoeffByOuterInner(int outer, int inner, const DenseBase<OtherDerived>& other);
+    template<typename OtherDerived>
     void copyCoeff(int index, const DenseBase<OtherDerived>& other);
     template<typename OtherDerived, int StoreMode, int LoadMode>
     void copyPacket(int row, int col, const DenseBase<OtherDerived>& other);
     template<typename OtherDerived, int StoreMode, int LoadMode>
+    void copyPacketByOuterInner(int outer, int inner, const DenseBase<OtherDerived>& other);
+    template<typename OtherDerived, int StoreMode, int LoadMode>
     void copyPacket(int index, const DenseBase<OtherDerived>& other);
+
+  private:
+    static int rowIndexByOuterInner(int outer, int inner);
+    static int colIndexByOuterInner(int outer, int inner);
+  public:
 #endif // not EIGEN_PARSED_BY_DOXYGEN
 
     template<int LoadMode>
     PacketScalar packet(int row, int col) const;
+    template<int LoadMode>
+    PacketScalar packetByOuterInner(int outer, int inner) const;
     template<int StoreMode>
     void writePacket(int row, int col, const PacketScalar& x);
+    template<int StoreMode>
+    void writePacketByOuterInner(int outer, int inner, const PacketScalar& x);
 
     template<int LoadMode>
     PacketScalar packet(int index) const;
@@ -409,13 +500,6 @@ template<typename Derived> class DenseBase
     template<typename OtherDerived>
     void swap(DenseBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other);
 
-    /** \returns number of elements to skip to pass from one row (resp. column) to another
-      * for a row-major (resp. column-major) matrix.
-      * Combined with coeffRef() and the \ref flags flags, it allows a direct access to the data
-      * of the underlying matrix.
-      */
-    inline int stride() const { return derived().stride(); }
-
     inline const NestByValue<Derived> nestByValue() const;
     inline const ForceAlignedAccess<Derived> forceAlignedAccess() const;
     inline ForceAlignedAccess<Derived> forceAlignedAccess();
diff --git a/Eigen/src/Core/DenseStorageBase.h b/Eigen/src/Core/DenseStorageBase.h
index 6245b8007..89e6e7112 100644
--- a/Eigen/src/Core/DenseStorageBase.h
+++ b/Eigen/src/Core/DenseStorageBase.h
@@ -75,23 +75,6 @@ class DenseStorageBase : public _Base<Derived>
     EIGEN_STRONG_INLINE int rows() const { return m_storage.rows(); }
     EIGEN_STRONG_INLINE int cols() const { return m_storage.cols(); }
 
-    /** 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 data() */
-    EIGEN_STRONG_INLINE int stride() const
-    {
-      if(IsVectorAtCompileTime)
-        return 1;
-      else
-        return (Flags & RowMajorBit) ? m_storage.cols() : m_storage.rows();
-    }
-
     EIGEN_STRONG_INLINE const Scalar& coeff(int row, int col) const
     {
       if(Flags & RowMajorBit)
@@ -253,12 +236,12 @@ class DenseStorageBase : public _Base<Derived>
     {
       if(RowsAtCompileTime == 1)
       {
-        ei_assert(other.isVector());
+        ei_assert(other.rows() == 1 || other.cols() == 1);
         resize(1, other.size());
       }
       else if(ColsAtCompileTime == 1)
       {
-        ei_assert(other.isVector());
+        ei_assert(other.rows() == 1 || other.cols() == 1);
         resize(other.size(), 1);
       }
       else resize(other.rows(), other.cols());
diff --git a/Eigen/src/Core/Flagged.h b/Eigen/src/Core/Flagged.h
index 0044fe7cb..9413b74fa 100644
--- a/Eigen/src/Core/Flagged.h
+++ b/Eigen/src/Core/Flagged.h
@@ -60,7 +60,8 @@ template<typename ExpressionType, unsigned int Added, unsigned int Removed> clas
 
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
-    inline int stride() const { return m_matrix.stride(); }
+    inline int outerStride() const { return m_matrix.outerStride(); }
+    inline int innerStride() const { return m_matrix.innerStride(); }
 
     inline const Scalar coeff(int row, int col) const
     {
diff --git a/Eigen/src/Core/ForceAlignedAccess.h b/Eigen/src/Core/ForceAlignedAccess.h
index 927f43413..300b22329 100644
--- a/Eigen/src/Core/ForceAlignedAccess.h
+++ b/Eigen/src/Core/ForceAlignedAccess.h
@@ -52,7 +52,8 @@ template<typename ExpressionType> class ForceAlignedAccess
 
     inline int rows() const { return m_expression.rows(); }
     inline int cols() const { return m_expression.cols(); }
-    inline int stride() const { return m_expression.stride(); }
+    inline int outerStride() const { return m_expression.outerStride(); }
+    inline int innerStride() const { return m_expression.innerStride(); }
 
     inline const CoeffReturnType coeff(int row, int col) const
     {
diff --git a/Eigen/src/Core/Map.h b/Eigen/src/Core/Map.h
index 83688dbca..d9ccb1b20 100644
--- a/Eigen/src/Core/Map.h
+++ b/Eigen/src/Core/Map.h
@@ -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) 2007-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
@@ -48,32 +48,53 @@
   *
   * \sa Matrix::Map()
   */
-template<typename MatrixType, int Options>
-struct ei_traits<Map<MatrixType, Options> > : public ei_traits<MatrixType>
+template<typename MatrixType, int Options, typename StrideType>
+struct ei_traits<Map<MatrixType, Options, StrideType> >
+  : public ei_traits<MatrixType>
 {
   enum {
-    Flags = (Options&Aligned)==Aligned ? ei_traits<MatrixType>::Flags |  AlignedBit
-                                       : ei_traits<MatrixType>::Flags & ~AlignedBit
+    Flags0 = ei_traits<MatrixType>::Flags,
+    Flags1 = ((Options&Aligned)==Aligned ? Flags0 |  AlignedBit
+                                         : Flags0 & ~AlignedBit),
+    Flags = int(StrideType::InnerStrideAtCompileTime)==1 ? Flags1 : (Flags1 & ~PacketAccessBit)
   };
 };
 
-template<typename MatrixType, int Options> class Map
-  : public MapBase<Map<MatrixType, Options>,
-                   typename MatrixType::template MakeBase< Map<MatrixType, Options> >::Type>
+template<typename MatrixType, int Options, typename StrideType> class Map
+  : public MapBase<Map<MatrixType, Options, StrideType>,
+                   typename MatrixType::template MakeBase<
+                     Map<MatrixType, Options, StrideType>
+                   >::Type>
 {
   public:
 
     typedef MapBase<Map,typename MatrixType::template MakeBase<Map>::Type> Base;
+
     EIGEN_DENSE_PUBLIC_INTERFACE(Map)
 
-    inline int stride() const { return this->innerSize(); }
+    inline int innerStride() const
+    {
+      return StrideType::InnerStrideAtCompileTime != 0 ? m_stride.inner() : 1;
+    }
 
-    inline Map(const Scalar* data) : Base(data) {}
+    inline int outerStride() const
+    {
+      return StrideType::OuterStrideAtCompileTime != 0 ? m_stride.outer()
+           : IsVectorAtCompileTime ? this->size()
+           : int(Flags)&RowMajorBit ? this->cols()
+           : this->rows();
+    }
+
+    inline Map(const Scalar* data, const StrideType& stride = StrideType())
+      : Base(data), m_stride(stride) {}
 
-    inline Map(const Scalar* data, int size) : Base(data, size) {}
+    inline Map(const Scalar* data, int size, const StrideType& stride = StrideType())
+      : Base(data, size), m_stride(stride) {}
 
-    inline Map(const Scalar* data, int rows, int cols) : Base(data, rows, cols) {}
+    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);
@@ -88,8 +109,12 @@ template<typename MatrixType, int Options> class Map
       EIGEN_ONLY_USED_FOR_DEBUG(size);
       ei_assert(size == this->size());
     }
+*/
 
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Map)
+
+  protected:
+    StrideType m_stride;
 };
 
 template<typename _Scalar, int _Rows, int _Cols, int _StorageOrder, int _MaxRows, int _MaxCols>
diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h
index a922d8bb0..d735cfc47 100644
--- a/Eigen/src/Core/MapBase.h
+++ b/Eigen/src/Core/MapBase.h
@@ -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) 2007-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
@@ -37,9 +37,7 @@ template<typename Derived, typename Base> class MapBase
 {
   public:
 
-//     typedef MatrixBase<Derived> Base;
     enum {
-      IsRowMajor = (int(ei_traits<Derived>::Flags) & RowMajorBit) ? 1 : 0,
       RowsAtCompileTime = ei_traits<Derived>::RowsAtCompileTime,
       ColsAtCompileTime = ei_traits<Derived>::ColsAtCompileTime,
       SizeAtCompileTime = Base::SizeAtCompileTime
@@ -48,94 +46,75 @@ template<typename Derived, typename Base> class MapBase
     typedef typename ei_traits<Derived>::Scalar Scalar;
     typedef typename Base::PacketScalar PacketScalar;
     using Base::derived;
+    using Base::innerStride;
+    using Base::outerStride;
+    using Base::rowStride;
+    using Base::colStride;
 
     inline int rows() const { return m_rows.value(); }
     inline int cols() const { return m_cols.value(); }
 
-    /** 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 MapBase::data() */
-    inline int stride() const { return derived().stride(); }
-
     /** Returns a pointer to the first coefficient of the matrix or vector.
-      * This function has to be used together with the stride() function.
       *
-      * \sa MapBase::stride() */
+      * \note When addressing this data, make sure to honor the strides returned by innerStride() and outerStride().
+      *
+      * \sa innerStride(), outerStride()
+      */
     inline const Scalar* data() const { return m_data; }
 
     inline const Scalar& coeff(int row, int col) const
     {
-      if(IsRowMajor)
-        return m_data[col + row * stride()];
-      else // column-major
-        return m_data[row + col * stride()];
+      return m_data[col * colStride() + row * rowStride()];
     }
 
     inline Scalar& coeffRef(int row, int col)
     {
-      if(IsRowMajor)
-        return const_cast<Scalar*>(m_data)[col + row * stride()];
-      else // column-major
-        return const_cast<Scalar*>(m_data)[row + col * stride()];
+      return const_cast<Scalar*>(m_data)[col * colStride() + row * rowStride()];
     }
 
     inline const Scalar& coeff(int index) const
     {
       ei_assert(Derived::IsVectorAtCompileTime || (ei_traits<Derived>::Flags & LinearAccessBit));
-      if ( ((RowsAtCompileTime == 1) == IsRowMajor) || !int(Derived::IsVectorAtCompileTime) )
-        return m_data[index];
-      else
-        return m_data[index*stride()];
+      return m_data[index * innerStride()];
     }
 
     inline Scalar& coeffRef(int index)
     {
       ei_assert(Derived::IsVectorAtCompileTime || (ei_traits<Derived>::Flags & LinearAccessBit));
-      if ( ((RowsAtCompileTime == 1) == IsRowMajor) || !int(Derived::IsVectorAtCompileTime) )
-        return const_cast<Scalar*>(m_data)[index];
-      else
-        return const_cast<Scalar*>(m_data)[index*stride()];
+      return const_cast<Scalar*>(m_data)[index * innerStride()];
     }
 
     template<int LoadMode>
     inline PacketScalar packet(int row, int col) const
     {
       return ei_ploadt<Scalar, LoadMode>
-               (m_data + (IsRowMajor ? col + row * stride()
-                                     : row + col * stride()));
+               (m_data + (col * colStride() + row * rowStride()));
     }
 
     template<int LoadMode>
     inline PacketScalar packet(int index) const
     {
-      return ei_ploadt<Scalar, LoadMode>(m_data + index);
+      return ei_ploadt<Scalar, LoadMode>(m_data + index * innerStride());
     }
 
     template<int StoreMode>
     inline void writePacket(int row, int col, const PacketScalar& x)
     {
       ei_pstoret<Scalar, PacketScalar, StoreMode>
-               (const_cast<Scalar*>(m_data) + (IsRowMajor ? col + row * stride()
-                                                          : row + col * stride()), x);
+               (const_cast<Scalar*>(m_data) + (col * colStride() + row * rowStride()), x);
     }
 
     template<int StoreMode>
     inline void writePacket(int index, const PacketScalar& x)
     {
       ei_pstoret<Scalar, PacketScalar, StoreMode>
-        (const_cast<Scalar*>(m_data) + index, x);
+        (const_cast<Scalar*>(m_data) + index * innerStride(), x);
     }
 
     inline MapBase(const Scalar* data) : m_data(data), m_rows(RowsAtCompileTime), m_cols(ColsAtCompileTime)
     {
       EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived)
-      checkDataAlignment();
+      checkSanity();
     }
 
     inline MapBase(const Scalar* data, int size)
@@ -146,7 +125,7 @@ template<typename Derived, typename Base> class MapBase
       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
       ei_assert(size >= 0);
       ei_assert(data == 0 || SizeAtCompileTime == Dynamic || SizeAtCompileTime == size);
-      checkDataAlignment();
+      checkSanity();
     }
 
     inline MapBase(const Scalar* data, int rows, int cols)
@@ -155,7 +134,7 @@ template<typename Derived, typename Base> class MapBase
       ei_assert( (data == 0)
               || (   rows >= 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
                   && cols >= 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols)));
-      checkDataAlignment();
+      checkSanity();
     }
 
     Derived& operator=(const MapBase& other)
@@ -167,10 +146,12 @@ template<typename Derived, typename Base> class MapBase
 
   protected:
 
-    void checkDataAlignment() const
+    void checkSanity() const
     {
       ei_assert( ((!(ei_traits<Derived>::Flags&AlignedBit))
                   || ((size_t(m_data)&0xf)==0)) && "data is not aligned");
+      ei_assert( ((!(ei_traits<Derived>::Flags&PacketAccessBit))
+                  || (innerStride()==1)) && "packet access incompatible with inner stride greater than 1");
     }
 
     const Scalar* EIGEN_RESTRICT m_data;
diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index 44a0ef7d1..eae2711f4 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -318,6 +318,9 @@ class Matrix
     void swap(MatrixBase<OtherDerived> EIGEN_REF_TO_TEMPORARY other)
     { this->_swap(other.derived()); }
 
+    inline int innerStride() const { return 1; }
+    inline int outerStride() const { return this->innerSize(); }
+
     /////////// Geometry module ///////////
 
     template<typename OtherDerived>
diff --git a/Eigen/src/Core/NestByValue.h b/Eigen/src/Core/NestByValue.h
index 9f6d1c0c0..497ce828b 100644
--- a/Eigen/src/Core/NestByValue.h
+++ b/Eigen/src/Core/NestByValue.h
@@ -53,7 +53,8 @@ template<typename ExpressionType> class NestByValue
 
     inline int rows() const { return m_expression.rows(); }
     inline int cols() const { return m_expression.cols(); }
-    inline int stride() const { return m_expression.stride(); }
+    inline int outerStride() const { return m_expression.outerStride(); }
+    inline int innerStride() const { return m_expression.innerStride(); }
 
     inline const CoeffReturnType coeff(int row, int col) const
     {
diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h
index d375f0b5c..160b973bd 100644
--- a/Eigen/src/Core/ReturnByValue.h
+++ b/Eigen/src/Core/ReturnByValue.h
@@ -34,14 +34,9 @@ struct ei_traits<ReturnByValue<Derived> >
   : public ei_traits<typename ei_traits<Derived>::ReturnType>
 {
   enum {
-    // FIXME had to remove the DirectAccessBit for usage like
-    //   matrix.inverse().block(...)
-    // because the Block ctor with direct access
-    // wants to call coeffRef() to get an address, and that fails (infinite recursion) as ReturnByValue
-    // doesnt implement coeffRef().
-    // The fact that I had to do that shows that when doing xpr.block() with a non-direct-access xpr,
-    // even if xpr has the EvalBeforeNestingBit, the block() doesn't use direct access on the evaluated
-    // xpr.
+    // We're disabling the DirectAccess because e.g. the constructor of
+    // the Block-with-DirectAccess expression requires to have a coeffRef method.
+    // Also, we don't want to have to implement the stride stuff.
     Flags = (ei_traits<typename ei_traits<Derived>::ReturnType>::Flags
              | EvalBeforeNestingBit) & ~DirectAccessBit
   };
diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
index c3c5b17ff..0b57b9968 100644
--- a/Eigen/src/Core/SelfAdjointView.h
+++ b/Eigen/src/Core/SelfAdjointView.h
@@ -75,8 +75,9 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
 
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
-    inline int stride() const { return m_matrix.stride(); }
-
+    inline int outerStride() const { return m_matrix.outerStride(); }
+    inline int innerStride() const { return m_matrix.innerStride(); }
+    
     /** \sa MatrixBase::coeff()
       * \warning the coordinates must fit into the referenced triangular part
       */
diff --git a/Eigen/src/Core/SelfCwiseBinaryOp.h b/Eigen/src/Core/SelfCwiseBinaryOp.h
index d2690b66b..529a9994d 100644
--- a/Eigen/src/Core/SelfCwiseBinaryOp.h
+++ b/Eigen/src/Core/SelfCwiseBinaryOp.h
@@ -57,7 +57,8 @@ template<typename BinaryOp, typename MatrixType> class SelfCwiseBinaryOp
 
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
-    inline int stride() const { return m_matrix.stride(); }
+    inline int outerStride() const { return m_matrix.outerStride(); }
+    inline int innerStride() const { return m_matrix.innerStride(); }
     inline const Scalar* data() const { return m_matrix.data(); }
 
     // note that this function is needed by assign to correctly align loads/stores
diff --git a/Eigen/src/Core/Stride.h b/Eigen/src/Core/Stride.h
new file mode 100644
index 000000000..7982035fd
--- /dev/null
+++ b/Eigen/src/Core/Stride.h
@@ -0,0 +1,80 @@
+// 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/>.
+
+#ifndef EIGEN_STRIDE_H
+#define EIGEN_STRIDE_H
+
+template<int _InnerStrideAtCompileTime, int _OuterStrideAtCompileTime>
+class Stride
+{
+  public:
+
+    enum {
+      InnerStrideAtCompileTime = _InnerStrideAtCompileTime,
+      OuterStrideAtCompileTime = _OuterStrideAtCompileTime
+    };
+
+    Stride()
+      : m_inner(InnerStrideAtCompileTime), m_outer(OuterStrideAtCompileTime)
+    {
+      ei_assert(InnerStrideAtCompileTime != Dynamic && OuterStrideAtCompileTime != Dynamic);
+    }
+
+    Stride(int innerStride, int outerStride)
+      : m_inner(innerStride), m_outer(outerStride)
+    {
+      ei_assert(innerStride>=0 && outerStride>=0);
+    }
+
+    Stride(const Stride& other)
+      : m_inner(other.inner()), m_outer(other.outer())
+    {}
+
+    inline int inner() const { return m_inner.value(); }
+    inline int outer() const { return m_outer.value(); }
+
+  protected:
+    ei_int_if_dynamic<InnerStrideAtCompileTime> m_inner;
+    ei_int_if_dynamic<OuterStrideAtCompileTime> m_outer;
+};
+
+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>
+{
+    typedef Stride<0,Value> Base;
+  public:
+    OuterStride() : Base() {}
+    OuterStride(int v) : Base(0,v) {}
+};
+
+#endif // EIGEN_STRIDE_H
diff --git a/Eigen/src/Core/Swap.h b/Eigen/src/Core/Swap.h
index 186268af0..c3c641097 100644
--- a/Eigen/src/Core/Swap.h
+++ b/Eigen/src/Core/Swap.h
@@ -47,7 +47,8 @@ template<typename ExpressionType> class SwapWrapper
 
     inline int rows() const { return m_expression.rows(); }
     inline int cols() const { return m_expression.cols(); }
-    inline int stride() const { return m_expression.stride(); }
+    inline int outerStride() const { return m_expression.outerStride(); }
+    inline int innerStride() const { return m_expression.innerStride(); }
 
     inline Scalar& coeffRef(int row, int col)
     {
@@ -60,7 +61,7 @@ template<typename ExpressionType> class SwapWrapper
     }
 
     template<typename OtherDerived>
-    void copyCoeff(int row, int col, const MatrixBase<OtherDerived>& other)
+    void copyCoeff(int row, int col, const DenseBase<OtherDerived>& other)
     {
       OtherDerived& _other = other.const_cast_derived();
       ei_internal_assert(row >= 0 && row < rows()
@@ -71,7 +72,7 @@ template<typename ExpressionType> class SwapWrapper
     }
 
     template<typename OtherDerived>
-    void copyCoeff(int index, const MatrixBase<OtherDerived>& other)
+    void copyCoeff(int index, const DenseBase<OtherDerived>& other)
     {
       OtherDerived& _other = other.const_cast_derived();
       ei_internal_assert(index >= 0 && index < m_expression.size());
@@ -81,7 +82,7 @@ template<typename ExpressionType> class SwapWrapper
     }
 
     template<typename OtherDerived, int StoreMode, int LoadMode>
-    void copyPacket(int row, int col, const MatrixBase<OtherDerived>& other)
+    void copyPacket(int row, int col, const DenseBase<OtherDerived>& other)
     {
       OtherDerived& _other = other.const_cast_derived();
       ei_internal_assert(row >= 0 && row < rows()
@@ -94,7 +95,7 @@ template<typename ExpressionType> class SwapWrapper
     }
 
     template<typename OtherDerived, int StoreMode, int LoadMode>
-    void copyPacket(int index, const MatrixBase<OtherDerived>& other)
+    void copyPacket(int index, const DenseBase<OtherDerived>& other)
     {
       OtherDerived& _other = other.const_cast_derived();
       ei_internal_assert(index >= 0 && index < m_expression.size());
diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h
index 6c0e50de2..1f064d1c2 100644
--- a/Eigen/src/Core/Transpose.h
+++ b/Eigen/src/Core/Transpose.h
@@ -93,7 +93,8 @@ template<typename MatrixType> class TransposeImpl<MatrixType,Dense>
     typedef typename MatrixType::template MakeBase<Transpose<MatrixType> >::Type Base;
     EIGEN_DENSE_PUBLIC_INTERFACE(Transpose<MatrixType>)
 
-    inline int stride() const { return derived().nestedExpression().stride(); }
+    inline int innerStride() const { return derived().nestedExpression().innerStride(); }
+    inline int outerStride() const { return derived().nestedExpression().outerStride(); }
     inline Scalar* data() { return derived().nestedExpression().data(); }
     inline const Scalar* data() const { return derived().nestedExpression().data(); }
 
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index 7d978b800..2230680d1 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -50,7 +50,8 @@ template<typename Derived> class TriangularBase : public EigenBase<Derived>
 
     inline int rows() const { return derived().rows(); }
     inline int cols() const { return derived().cols(); }
-    inline int stride() const { return derived().stride(); }
+    inline int outerStride() const { return derived().outerStride(); }
+    inline int innerStride() const { return derived().innerStride(); }
 
     inline Scalar coeff(int row, int col) const  { return derived().coeff(row,col); }
     inline Scalar& coeffRef(int row, int col) { return derived().coeffRef(row,col); }
@@ -165,7 +166,8 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
 
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
-    inline int stride() const { return m_matrix.stride(); }
+    inline int outerStride() const { return m_matrix.outerStride(); }
+    inline int innerStride() const { return m_matrix.innerStride(); }
 
     /** \sa MatrixBase::operator+=() */
     template<typename Other> TriangularView&  operator+=(const Other& other) { return *this = m_matrix + other; }
diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h
index cbf97aeb3..5bb7fd35d 100644
--- a/Eigen/src/Core/VectorBlock.h
+++ b/Eigen/src/Core/VectorBlock.h
@@ -86,7 +86,6 @@ template<typename VectorType, int Size> class VectorBlock
              IsColVector ? start : 0, IsColVector ? 0 : start,
              IsColVector ? size  : 1, IsColVector ? 1 : size)
     {
-
       EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock);
     }
 
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 51590b03d..c27c979a6 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -140,7 +140,7 @@ const unsigned int LinearAccessBit = 0x10;
   * Means that the underlying array of coefficients can be directly accessed. This means two things.
   * First, references to the coefficients must be available through coeffRef(int, int). This rules out read-only
   * expressions whose coefficients are computed on demand by coeff(int, int). Second, the memory layout of the
-  * array of coefficients must be exactly the natural one suggested by rows(), cols(), stride(), and the RowMajorBit.
+  * array of coefficients must be exactly the natural one suggested by rows(), cols(), outerStride(), innerStride(), and the RowMajorBit.
   * This rules out expressions such as Diagonal, whose coefficients, though referencable, do not have
   * such a regular memory layout.
   */
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index a56e4fb84..6096272fa 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -1,7 +1,8 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2007-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -60,7 +61,9 @@ 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<typename MatrixType, int Options=Unaligned> class Map;
+template<int InnerStrideAtCompileTime = Dynamic, int OuterStrideAtCompileTime = Dynamic> class Stride;
+template<typename MatrixType, int Options=Unaligned, typename StrideType = Stride<0,0> > class Map;
+
 template<typename Derived> class TriangularBase;
 template<typename MatrixType, unsigned int Mode> class TriangularView;
 template<typename MatrixType, unsigned int Mode> class SelfAdjointView;
diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h
index d4920d213..c7b95d334 100644
--- a/Eigen/src/Core/util/Memory.h
+++ b/Eigen/src/Core/util/Memory.h
@@ -232,7 +232,7 @@ inline static Integer ei_first_aligned(const Scalar* array, Integer size)
   enum { PacketSize = ei_packet_traits<Scalar>::size,
          PacketAlignedMask = PacketSize-1
   };
-  
+
   if(PacketSize==1)
   {
     // Either there is no vectorization, or a packet consists of exactly 1 scalar so that all elements
diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h
index 619e7664d..5252b28c5 100644
--- a/Eigen/src/Core/util/StaticAssert.h
+++ b/Eigen/src/Core/util/StaticAssert.h
@@ -2,7 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2008 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
@@ -46,7 +46,7 @@
     // if native static_assert is enabled, let's use it
     #define EIGEN_STATIC_ASSERT(X,MSG) static_assert(X,#MSG);
 
-  #else // CXX0X
+  #else // not CXX0X
 
     template<bool condition>
     struct ei_static_assert {};
@@ -81,7 +81,8 @@
         BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER,
         THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX,
         THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE,
-        THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES
+        THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES,
+        YOU_ALREADY_SPECIFIED_THIS_STRIDE
       };
     };
 
diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h
index f1ed4ed9f..eceb5ab2a 100644
--- a/Eigen/src/Core/util/XprHelper.h
+++ b/Eigen/src/Core/util/XprHelper.h
@@ -50,7 +50,7 @@ template<int Value> class ei_int_if_dynamic
 {
   public:
     EIGEN_EMPTY_STRUCT_CTOR(ei_int_if_dynamic)
-    explicit ei_int_if_dynamic(int) {}
+    explicit ei_int_if_dynamic(int v) { EIGEN_ONLY_USED_FOR_DEBUG(v); ei_assert(v == Value); }
     static int value() { return Value; }
     void setValue(int) {}
 };
@@ -58,7 +58,7 @@ template<int Value> class ei_int_if_dynamic
 template<> class ei_int_if_dynamic<Dynamic>
 {
     int m_value;
-    ei_int_if_dynamic() {}
+    ei_int_if_dynamic() { ei_assert(false); }
   public:
     explicit ei_int_if_dynamic(int value) : m_value(value) {}
     int value() const { return m_value; }