* Add fixed-size template versions of corner(), start(), end().

* Use them to write an unrolled path in echelon.cpp, as an
  experiment before I do this LU module.
* For floating-point types, make ei_random() use an amplitude
  of 1.

Eigen/src/Core/Block.h

@@ -67,9 +67,11 @@ struct ei_traits<Block<MatrixType, BlockRows, BlockCols> >
       : (BlockRows==Dynamic ? MatrixType::MaxRowsAtCompileTime : BlockRows),
     MaxColsAtCompileTime = ColsAtCompileTime == 1 ? 1
       : (BlockCols==Dynamic ? MatrixType::MaxColsAtCompileTime : BlockCols),
-    Flags = (RowsAtCompileTime == Dynamic || ColsAtCompileTime == Dynamic
-            ? (unsigned int)MatrixType::Flags
-            : (unsigned int)MatrixType::Flags &~ LargeBit) & ~VectorizableBit,
+    FlagsLargeBit = ((RowsAtCompileTime != Dynamic && MatrixType::RowsAtCompileTime == Dynamic)
+                  || (ColsAtCompileTime != Dynamic && MatrixType::ColsAtCompileTime == Dynamic))
+                  ? ~LargeBit
+                  : ~(unsigned int)0,
+    Flags = MatrixType::Flags & FlagsLargeBit & ~VectorizableBit,
     CoeffReadCost = MatrixType::CoeffReadCost
   };
 };
@@ -236,8 +238,7 @@ const Block<Derived> MatrixBase<Derived>
   * \sa class Block, block(int,int)
   */
 template<typename Derived>
-Block<Derived> MatrixBase<Derived>
-  ::start(int size)
+Block<Derived> MatrixBase<Derived>::start(int size)
 {
   ei_assert(IsVectorAtCompileTime);
   return Block<Derived>(derived(), 0, 0,
@@ -247,8 +248,7 @@ Block<Derived> MatrixBase<Derived>
 
 /** This is the const version of start(int).*/
 template<typename Derived>
-const Block<Derived> MatrixBase<Derived>
-  ::start(int size) const
+const Block<Derived> MatrixBase<Derived>::start(int size) const
 {
   ei_assert(IsVectorAtCompileTime);
   return Block<Derived>(derived(), 0, 0,
@@ -272,8 +272,7 @@ const Block<Derived> MatrixBase<Derived>
   * \sa class Block, block(int,int)
   */
 template<typename Derived>
-Block<Derived> MatrixBase<Derived>
-  ::end(int size)
+Block<Derived> MatrixBase<Derived>::end(int size)
 {
   ei_assert(IsVectorAtCompileTime);
   return Block<Derived>(derived(),
@@ -285,8 +284,7 @@ Block<Derived> MatrixBase<Derived>
 
 /** This is the const version of end(int).*/
 template<typename Derived>
-const Block<Derived> MatrixBase<Derived>
-  ::end(int size) const
+const Block<Derived> MatrixBase<Derived>::end(int size) const
 {
   ei_assert(IsVectorAtCompileTime);
   return Block<Derived>(derived(),
@@ -296,6 +294,80 @@ const Block<Derived> MatrixBase<Derived>
                         ColsAtCompileTime == 1 ? 1 : size);
 }
 
+/** \returns a fixed-size expression of the first coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * The template parameter \a Size is the number of coefficients in the block
+  *
+  * Example: \include MatrixBase_template_int_start.cpp
+  * Output: \verbinclude MatrixBase_template_int_start.out
+  *
+  * \sa class Block
+  */
+template<typename Derived>
+template<int Size>
+Block<Derived, ei_traits<Derived>::RowsAtCompileTime == 1 ? 1 : Size,
+               ei_traits<Derived>::ColsAtCompileTime == 1 ? 1 : Size>
+MatrixBase<Derived>::start()
+{
+  ei_assert(IsVectorAtCompileTime);
+  return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size,
+                        ColsAtCompileTime == 1 ? 1 : Size>(derived(), 0, 0);
+}
+
+/** This is the const version of start<int>().*/
+template<typename Derived>
+template<int Size>
+const Block<Derived, ei_traits<Derived>::RowsAtCompileTime == 1 ? 1 : Size,
+                     ei_traits<Derived>::ColsAtCompileTime == 1 ? 1 : Size>
+MatrixBase<Derived>::start() const
+{
+  ei_assert(IsVectorAtCompileTime);
+  return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size,
+                        ColsAtCompileTime == 1 ? 1 : Size>(derived(), 0, 0);
+}
+
+/** \returns a fixed-size expression of the last coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * The template parameter \a Size is the number of coefficients in the block
+  *
+  * Example: \include MatrixBase_template_int_end.cpp
+  * Output: \verbinclude MatrixBase_template_int_end.out
+  *
+  * \sa class Block
+  */
+template<typename Derived>
+template<int Size>
+Block<Derived, ei_traits<Derived>::RowsAtCompileTime == 1 ? 1 : Size,
+               ei_traits<Derived>::ColsAtCompileTime == 1 ? 1 : Size>
+MatrixBase<Derived>::end()
+{
+  ei_assert(IsVectorAtCompileTime);
+  return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size,
+                        ColsAtCompileTime == 1 ? 1 : Size>
+           (derived(),
+            RowsAtCompileTime == 1 ? 0 : rows() - Size,
+            ColsAtCompileTime == 1 ? 0 : cols() - Size);
+}
+
+/** This is the const version of end<int>.*/
+template<typename Derived>
+template<int Size>
+const Block<Derived, ei_traits<Derived>::RowsAtCompileTime == 1 ? 1 : Size,
+                     ei_traits<Derived>::ColsAtCompileTime == 1 ? 1 : Size>
+MatrixBase<Derived>::end() const
+{
+  ei_assert(IsVectorAtCompileTime);
+  return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size,
+                        ColsAtCompileTime == 1 ? 1 : Size>
+           (derived(),
+            RowsAtCompileTime == 1 ? 0 : rows() - Size,
+            ColsAtCompileTime == 1 ? 0 : cols() - Size);
+}
+
 /** \returns a dynamic-size expression of a corner of *this.
   *
   * \param type the type of corner. Can be \a Eigen::TopLeft, \a Eigen::TopRight,
@@ -316,14 +388,23 @@ template<typename Derived>
 Block<Derived> MatrixBase<Derived>
   ::corner(CornerType type, int cRows, int cCols)
 {
-  if(type == TopLeft)
-    return Block<Derived>(derived(), 0, 0, cRows, cCols);
-  else if(type == TopRight)
-    return Block<Derived>(derived(), 0, cols() - cCols, cRows, cCols);
-  else if(type == BottomLeft)
-    return Block<Derived>(derived(), rows() - cRows, 0, cRows, cCols);
-  else
-    return Block<Derived>(derived(), rows() - cRows, cols() - cCols, cRows, cCols);
+  switch(type)
+  {
+    case TopLeft:
+      return Block<Derived>(derived(), 0, 0, cRows, cCols);
+      break;
+    case TopRight:
+      return Block<Derived>(derived(), 0, cols() - cCols, cRows, cCols);
+      break;
+    case BottomLeft:
+      return Block<Derived>(derived(), rows() - cRows, 0, cRows, cCols);
+      break;
+    case BottomRight:
+      return Block<Derived>(derived(), rows() - cRows, cols() - cCols, cRows, cCols);
+      break;
+    default:
+      ei_assert(false && "Bad corner type.");
+  }
 }
 
 /** This is the const version of corner(CornerType, int, int).*/
@@ -331,14 +412,84 @@ template<typename Derived>
 const Block<Derived> MatrixBase<Derived>
   ::corner(CornerType type, int cRows, int cCols) const
 {
-  if(type == TopLeft)
-    return Block<Derived>(derived(), 0, 0, cRows, cCols);
-  else if(type == TopRight)
-    return Block<Derived>(derived(), 0, cols() - cCols, cRows, cCols);
-  else if(type == BottomLeft)
-    return Block<Derived>(derived(), rows() - cRows, 0, cRows, cCols);
-  else
-    return Block<Derived>(derived(), rows() - cRows, cols() - cCols, cRows, cCols);
+  switch(type)
+  {
+    case TopLeft:
+      return Block<Derived>(derived(), 0, 0, cRows, cCols);
+      break;
+    case TopRight:
+      return Block<Derived>(derived(), 0, cols() - cCols, cRows, cCols);
+      break;
+    case BottomLeft:
+      return Block<Derived>(derived(), rows() - cRows, 0, cRows, cCols);
+      break;
+    case BottomRight:
+      return Block<Derived>(derived(), rows() - cRows, cols() - cCols, cRows, cCols);
+      break;
+    default:
+      ei_assert(false && "Bad corner type.");
+  }
+}
+
+/** \returns a fixed-size expression of a corner of *this.
+  *
+  * \param type the type of corner. Can be \a Eigen::TopLeft, \a Eigen::TopRight,
+  * \a Eigen::BottomLeft, \a Eigen::BottomRight.
+  *
+  * The template parameters CRows and CCols arethe number of rows and columns in the corner.
+  *
+  * Example: \include MatrixBase_template_int_int_corner_enum.cpp
+  * Output: \verbinclude MatrixBase_template_int_int_corner_enum.out
+  *
+  * \sa class Block, block(int,int,int,int)
+  */
+template<typename Derived>
+template<int CRows, int CCols>
+Block<Derived, CRows, CCols> MatrixBase<Derived>
+  ::corner(CornerType type)
+{
+  switch(type)
+  {
+    case TopLeft:
+      return Block<Derived, CRows, CCols>(derived(), 0, 0);
+      break;
+    case TopRight:
+      return Block<Derived, CRows, CCols>(derived(), 0, cols() - CCols);
+      break;
+    case BottomLeft:
+      return Block<Derived, CRows, CCols>(derived(), rows() - CRows, 0);
+      break;
+    case BottomRight:
+      return Block<Derived, CRows, CCols>(derived(), rows() - CRows, cols() - CCols);
+      break;
+    default:
+      ei_assert(false && "Bad corner type.");
+  }
+}
+
+/** This is the const version of corner<int, int>(CornerType).*/
+template<typename Derived>
+template<int CRows, int CCols>
+const Block<Derived, CRows, CCols> MatrixBase<Derived>
+  ::corner(CornerType type) const
+{
+  switch(type)
+  {
+    case TopLeft:
+      return Block<Derived, CRows, CCols>(derived(), 0, 0);
+      break;
+    case TopRight:
+      return Block<Derived, CRows, CCols>(derived(), 0, cols() - CCols);
+      break;
+    case BottomLeft:
+      return Block<Derived, CRows, CCols>(derived(), rows() - CRows, 0);
+      break;
+    case BottomRight:
+      return Block<Derived, CRows, CCols>(derived(), rows() - CRows, cols() - CCols);
+      break;
+    default:
+      ei_assert(false && "Bad corner type.");
+  }
 }
 
 /** \returns a fixed-size expression of a block in *this.

Eigen/src/Core/MathFunctions.h

@@ -28,6 +28,11 @@
 template<typename T> inline typename NumTraits<T>::Real precision();
 template<typename T> inline T ei_random(T a, T b);
 template<typename T> inline T ei_random();
+template<typename T> inline T ei_random_amplitude()
+{
+  if(NumTraits<T>::HasFloatingPoint) return static_cast<T>(1);
+  else return static_cast<T>(10);
+}
 
 template<> inline int precision<int>() { return 0; }
 inline int ei_real(int x)  { return x; }
@@ -54,7 +59,7 @@ template<> inline int ei_random(int a, int b)
 }
 template<> inline int ei_random()
 {
-  return ei_random<int>(-10, 10);
+  return ei_random<int>(-ei_random_amplitude<int>(), ei_random_amplitude<int>());
 }
 inline bool ei_isMuchSmallerThan(int a, int, int = precision<int>())
 {
@@ -88,7 +93,7 @@ template<> inline float ei_random(float a, float b)
 }
 template<> inline float ei_random()
 {
-  return ei_random<float>(-10.0f, 10.0f);
+  return ei_random<float>(-ei_random_amplitude<float>(), ei_random_amplitude<float>());
 }
 inline bool ei_isMuchSmallerThan(float a, float b, float prec = precision<float>())
 {
@@ -122,7 +127,7 @@ template<> inline double ei_random(double a, double b)
 }
 template<> inline double ei_random()
 {
-  return ei_random<double>(-10.0, 10.0);
+  return ei_random<double>(-ei_random_amplitude<double>(), ei_random_amplitude<double>());
 }
 inline bool ei_isMuchSmallerThan(double a, double b, double prec = precision<double>())
 {

Eigen/src/Core/MatrixBase.h

@@ -316,6 +316,23 @@ template<typename Derived> class MatrixBase
     template<int BlockRows, int BlockCols>
     const Block<Derived, BlockRows, BlockCols> block(int startRow, int startCol) const;
 
+    template<int CRows, int CCols> Block<Derived, CRows, CCols> corner(CornerType type);
+    template<int CRows, int CCols> const Block<Derived, CRows, CCols> corner(CornerType type) const;
+
+    template<int Size>
+    Block<Derived, ei_traits<Derived>::RowsAtCompileTime == 1 ? 1 : Size,
+                   ei_traits<Derived>::ColsAtCompileTime == 1 ? 1 : Size> start();
+    template<int Size>
+    const Block<Derived, ei_traits<Derived>::RowsAtCompileTime == 1 ? 1 : Size,
+                         ei_traits<Derived>::ColsAtCompileTime == 1 ? 1 : Size> start() const;
+
+    template<int Size>
+    Block<Derived, ei_traits<Derived>::RowsAtCompileTime == 1 ? 1 : Size,
+                   ei_traits<Derived>::ColsAtCompileTime == 1 ? 1 : Size> end();
+    template<int Size>
+    const Block<Derived, ei_traits<Derived>::RowsAtCompileTime == 1 ? 1 : Size,
+                         ei_traits<Derived>::ColsAtCompileTime == 1 ? 1 : Size> end() const;
+
     DiagonalCoeffs<Derived> diagonal();
     const DiagonalCoeffs<Derived> diagonal() const;
     //@}

Eigen/src/Core/PacketMath.h

@@ -79,7 +79,7 @@ inline float  ei_pfirst(const __m128&  a) { return _mm_cvtss_f32(a); }
 inline double ei_pfirst(const __m128d& a) { return _mm_cvtsd_f64(a); }
 inline int    ei_pfirst(const __m128i& a) { return _mm_cvtsi128_si32(a); }
 
-#endif
+#endif // EIGEN_VECTORIZE_SSE
 
 #endif // EIGEN_PACKET_MATH_H
 

Eigen/src/Core/util/Macros.h

@@ -142,6 +142,7 @@ EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=)
 #define _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \
 typedef BaseClass Base; \
 typedef typename Eigen::ei_traits<Derived>::Scalar Scalar; \
+typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \
 typedef typename Base::PacketScalar PacketScalar; \
 typedef typename Eigen::ei_nested<Derived>::type Nested; \
 typedef typename Eigen::ei_eval<Derived>::type Eval; \

bench/benchmark.cpp

@@ -28,7 +28,7 @@ int main(int argc, char *argv[])
     asm("#begin");
     for(int a = 0; a < REPEAT; a++)
     {
-        m = I + 0.00005 * (m + m*m);
+        m = Matrix<SCALAR,MATSIZE,MATSIZE>::ones() + 0.00005 * (m + m*m);
     }
     asm("#end");
     cout << m << endl;

bench/benchmarkXcwise.cpp

@@ -10,24 +10,24 @@ USING_PART_OF_NAMESPACE_EIGEN
 #endif
 
 #ifndef MATSIZE
-#define MATSIZE 400
+#define MATSIZE 1000000
 #endif
 
 #ifndef REPEAT
-#define REPEAT 10000
+#define REPEAT 1000
 #endif
 
 int main(int argc, char *argv[])
 {
-	MATTYPE I = MATTYPE::ones(MATSIZE,MATSIZE);
-	MATTYPE m(MATSIZE,MATSIZE);
-	for(int i = 0; i < MATSIZE; i++) for(int j = 0; j < MATSIZE; j++)
+	MATTYPE I = MATTYPE::ones(MATSIZE,1);
+	MATTYPE m(MATSIZE,1);
+	for(int i = 0; i < MATSIZE; i++) for(int j = 0; j < 1; j++)
 	{
-		m(i,j) = 0.1 * (i+j+1)/(MATSIZE*MATSIZE);
+		m(i,j) = 0.1 * (i+j+1)/MATSIZE/MATSIZE;
 	}
 	for(int a = 0; a < REPEAT; a++)
 	{
-		m = I + 0.00005 * (m + m/4);
+		m = MATTYPE::ones(MATSIZE,1) + 0.00005 * (m.cwiseProduct(m) + m/4);
 	}
 	cout << m(0,0) << endl;
 	return 0;

doc/echelon.cpp

@@ -4,21 +4,72 @@ USING_PART_OF_NAMESPACE_EIGEN
 
 namespace Eigen {
 
-template<typename Derived>
-void echelon(MatrixBase<Derived>& m)
+/* Echelon a matrix in-place:
+ *
+ * Meta-Unrolled version, for small fixed-size matrices
+ */
+template<typename Derived, int Step>
+struct unroll_echelon
 {
-  for(int k = 0; k < m.diagonal().size(); k++)
+  enum { k = Step - 1,
+         Rows = Derived::RowsAtCompileTime,
+         Cols = Derived::ColsAtCompileTime,
+         CornerRows = Rows - k,
+         CornerCols = Cols - k
+  };
+  static void run(MatrixBase<Derived>& m)
   {
+    unroll_echelon<Derived, Step-1>::run(m);
     int rowOfBiggest, colOfBiggest;
-    int cornerRows = m.rows()-k, cornerCols = m.cols()-k;
-    m.corner(BottomRight, cornerRows, cornerCols)
+    m.template corner<CornerRows, CornerCols>(BottomRight)
      .cwiseAbs()
      .maxCoeff(&rowOfBiggest, &colOfBiggest);
     m.row(k).swap(m.row(k+rowOfBiggest));
     m.col(k).swap(m.col(k+colOfBiggest));
-    m.corner(BottomRight, cornerRows-1, cornerCols)
-      -= m.col(k).end(cornerRows-1) * (m.row(k).end(cornerCols) / m(k,k));
+    m.template corner<CornerRows-1, CornerCols>(BottomRight)
+      -= m.col(k).template end<CornerRows-1>()
+       * (m.row(k).template end<CornerCols>() / m(k,k));
   }
+};
+
+template<typename Derived>
+struct unroll_echelon<Derived, 0>
+{
+  static void run(MatrixBase<Derived>& m) {}
+};
+
+/* Echelon a matrix in-place:
+ *
+ * Non-unrolled version, for dynamic-size matrices.
+ * (this version works for all matrices, but in the fixed-size case the other
+ * version is faster).
+ */
+template<typename Derived>
+struct unroll_echelon<Derived, Dynamic>
+{
+  static void run(MatrixBase<Derived>& m)
+  {
+    for(int k = 0; k < m.diagonal().size(); k++)
+    {
+      int rowOfBiggest, colOfBiggest;
+      int cornerRows = m.rows()-k, cornerCols = m.cols()-k;
+      m.corner(BottomRight, cornerRows, cornerCols)
+      .cwiseAbs()
+      .maxCoeff(&rowOfBiggest, &colOfBiggest);
+      m.row(k).swap(m.row(k+rowOfBiggest));
+      m.col(k).swap(m.col(k+colOfBiggest));
+      m.corner(BottomRight, cornerRows-1, cornerCols)
+        -= m.col(k).end(cornerRows-1) * (m.row(k).end(cornerCols) / m(k,k));
+    }
+  }
+};
+using namespace std;
+template<typename Derived>
+void echelon(MatrixBase<Derived>& m)
+{
+  const int size = DiagonalCoeffs<Derived>::SizeAtCompileTime;
+  const bool unroll = size <= 4;
+  unroll_echelon<Derived, unroll ? size : Dynamic>::run(m);
 }
 
 template<typename Derived>
@@ -49,7 +100,7 @@ using namespace std;
 int main(int, char **)
 {
   srand((unsigned int)time(0));
-  const int Rows = 6, Cols = 4;
+  const int Rows = 6, Cols = 6;
   typedef Matrix<double, Rows, Cols> Mat;
   const int N = Rows < Cols ? Rows : Cols;