From abc8c010807f0706b80bc0ef13303b6485473a57 Mon Sep 17 00:00:00 2001
From: Hauke Heibel <hauke.heibel@gmail.com>
Date: Sat, 20 Feb 2010 15:53:57 +0100
Subject: Renamed PlainMatrixType to PlainObject (Array != Matrix). Renamed
 ReturnByValue::ReturnMatrixType ReturnByValue::ReturnType (again, Array !=
 Matrix).

---
 Eigen/src/LU/FullPivLU.h    | 6 +++---
 Eigen/src/LU/Inverse.h      | 4 ++--
 Eigen/src/LU/PartialPivLU.h | 8 ++++----
 3 files changed, 9 insertions(+), 9 deletions(-)

(limited to 'Eigen/src/LU')

diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 72e878223..1129293d5 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -630,7 +630,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
       return;
     }
 
-    typename Rhs::PlainMatrixType c(rhs().rows(), rhs().cols());
+    typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
 
     // Step 1
     c = dec().permutationP() * rhs();
@@ -670,10 +670,10 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
   * \sa class FullPivLU
   */
 template<typename Derived>
-inline const FullPivLU<typename MatrixBase<Derived>::PlainMatrixType>
+inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
 MatrixBase<Derived>::fullPivLu() const
 {
-  return FullPivLU<PlainMatrixType>(eval());
+  return FullPivLU<PlainObject>(eval());
 }
 
 #endif // EIGEN_LU_H
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index 36392c8d8..e20da70d6 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -238,7 +238,7 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 4>
 template<typename MatrixType>
 struct ei_traits<ei_inverse_impl<MatrixType> >
 {
-  typedef typename MatrixType::PlainMatrixType ReturnMatrixType;
+  typedef typename MatrixType::PlainObject ReturnType;
 };
 
 template<typename MatrixType>
@@ -327,7 +327,7 @@ inline void MatrixBase<Derived>::computeInverseAndDetWithCheck(
   typedef typename ei_meta_if<
     RowsAtCompileTime == 2,
     typename ei_cleantype<typename ei_nested<Derived, 2>::type>::type,
-    PlainMatrixType
+    PlainObject
   >::ret MatrixType;
   ei_compute_inverse_and_det_with_check<MatrixType, ResultType>::run
     (derived(), absDeterminantThreshold, inverse, determinant, invertible);
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index 809e4aad6..ed2354d78 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -442,10 +442,10 @@ struct ei_solve_retval<PartialPivLU<_MatrixType>, Rhs>
   * \sa class PartialPivLU
   */
 template<typename Derived>
-inline const PartialPivLU<typename MatrixBase<Derived>::PlainMatrixType>
+inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
 MatrixBase<Derived>::partialPivLu() const
 {
-  return PartialPivLU<PlainMatrixType>(eval());
+  return PartialPivLU<PlainObject>(eval());
 }
 
 /** \lu_module
@@ -457,10 +457,10 @@ MatrixBase<Derived>::partialPivLu() const
   * \sa class PartialPivLU
   */
 template<typename Derived>
-inline const PartialPivLU<typename MatrixBase<Derived>::PlainMatrixType>
+inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
 MatrixBase<Derived>::lu() const
 {
-  return PartialPivLU<PlainMatrixType>(eval());
+  return PartialPivLU<PlainObject>(eval());
 }
 
 #endif // EIGEN_PARTIALLU_H
-- 
cgit v1.2.3


From 71fecd23713ed728a5b94fd066b22fc92d122b9d Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Sat, 20 Feb 2010 18:19:34 +0100
Subject: add missing return value

---
 Eigen/src/LU/FullPivLU.h | 1 +
 1 file changed, 1 insertion(+)

(limited to 'Eigen/src/LU')

diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 1129293d5..9afc448cc 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -251,6 +251,7 @@ template<typename _MatrixType> class FullPivLU
     {
       m_usePrescribedThreshold = true;
       m_prescribedThreshold = threshold;
+      return *this;
     }
 
     /** Allows to come back to the default behavior, letting Eigen use its default formula for
-- 
cgit v1.2.3


From 7dc75380c101b9b4f3882f78fe6a5e9ae8963cac Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Tue, 23 Feb 2010 09:04:59 -0500
Subject: * FullPivLU: replace "remaining==0" termination condition (from
 Golub) by a fuzzy compare   (fixes lu test failures when testing solve()) *
 LU test: set appropriate threshold and limit the number of times that a
 specially tricky test   is run. (fixes lu test failures when testing rank()).
 * Tests: rename createRandomMatrixOfRank to createRandomProjectionOfRank

---
 Eigen/src/LU/FullPivLU.h |  6 +++++-
 test/inverse.cpp         |  2 +-
 test/lu.cpp              | 31 ++++++++++++++++++++++++++-----
 test/main.h              |  4 ++--
 test/qr_colpivoting.cpp  |  4 ++--
 test/qr_fullpivoting.cpp |  2 +-
 6 files changed, 37 insertions(+), 12 deletions(-)

(limited to 'Eigen/src/LU')

diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 9afc448cc..ec551645b 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -404,6 +404,7 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
 
   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
   m_maxpivot = RealScalar(0);
+  RealScalar cutoff(0);
 
   for(int k = 0; k < size; ++k)
   {
@@ -418,8 +419,11 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
     row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
     col_of_biggest_in_corner += k; // need to add k to them.
 
+    // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix
+    if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();
+
     // if the pivot (hence the corner) is exactly zero, terminate to avoid generating nan/inf values
-    if(biggest_in_corner == RealScalar(0))
+    if(ei_abs(biggest_in_corner) < cutoff)
     {
       // before exiting, make sure to initialize the still uninitialized transpositions
       // in a sane state without destroying what we already have.
diff --git a/test/inverse.cpp b/test/inverse.cpp
index 713caf4a6..3f6138e0c 100644
--- a/test/inverse.cpp
+++ b/test/inverse.cpp
@@ -42,7 +42,7 @@ template<typename MatrixType> void inverse(const MatrixType& m)
              m2(rows, cols),
              mzero = MatrixType::Zero(rows, cols),
              identity = MatrixType::Identity(rows, rows);
-  createRandomMatrixOfRank(rows,rows,rows,m1);
+  createRandomProjectionOfRank(rows,rows,rows,m1);
   m2 = m1.inverse();
   VERIFY_IS_APPROX(m1, m2.inverse() );
 
diff --git a/test/lu.cpp b/test/lu.cpp
index 568db8230..02f6ec805 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -28,7 +28,11 @@ using namespace std;
 
 template<typename MatrixType> void lu_non_invertible()
 {
+  static int times_called = 0;
+  times_called++;
+  
   typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
   /* this test covers the following files:
      LU.h
   */
@@ -64,9 +68,15 @@ template<typename MatrixType> void lu_non_invertible()
   
   MatrixType m1(rows, cols), m3(rows, cols2);
   CMatrixType m2(cols, cols2);
-  createRandomMatrixOfRank(rank, rows, cols, m1);
+  createRandomProjectionOfRank(rank, rows, cols, m1);
+
+  FullPivLU<MatrixType> lu;
+
+  // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank of projections.
+  // So it's not clear at all the epsilon should play any role there.
+  lu.setThreshold(RealScalar(0.01));
+  lu.compute(m1);
 
-  FullPivLU<MatrixType> lu(m1);
   // FIXME need better way to construct trapezoid matrices. extend triangularView to support rectangular.
   DynamicMatrixType u(rows,cols);
   for(int i = 0; i < rows; i++)
@@ -91,9 +101,20 @@ template<typename MatrixType> void lu_non_invertible()
   VERIFY(!lu.isSurjective());
   VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
   VERIFY(m1image.fullPivLu().rank() == rank);
-  DynamicMatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
-  sidebyside << m1, m1image;
-  VERIFY(sidebyside.fullPivLu().rank() == rank);
+
+  // The following test is damn hard to get to succeed over a large number of repetitions.
+  // We're checking that the image is indeed the image, i.e. adding it as new columns doesn't increase the rank.
+  // Since we've already tested rank() above, the point here is not to test rank(), it is to test image().
+  // Since image() is implemented in a very simple way that doesn't leave much room for choice, the occasional
+  // errors that we get here (one in 1e+4 repetitions roughly) are probably just a sign that it's a really
+  // hard test, so we just limit how many times it's run.
+  if(times_called < 100)
+  {
+    DynamicMatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
+    sidebyside << m1, m1image;
+    VERIFY(sidebyside.fullPivLu().rank() == rank);
+  }
+  
   m2 = CMatrixType::Random(cols,cols2);
   m3 = m1*m2;
   m2 = CMatrixType::Random(cols,cols2);
diff --git a/test/main.h b/test/main.h
index 64f70b394..6d296b2e3 100644
--- a/test/main.h
+++ b/test/main.h
@@ -148,7 +148,7 @@ namespace Eigen
 
 #define EIGEN_INTERNAL_DEBUGGING
 #define EIGEN_NICE_RANDOM
-#include <Eigen/QR> // required for createRandomMatrixOfRank
+#include <Eigen/QR> // required for createRandomProjectionOfRank
 
 
 #define VERIFY(a) do { if (!(a)) { \
@@ -343,7 +343,7 @@ inline bool test_isUnitary(const MatrixBase<Derived>& m)
 }
 
 template<typename MatrixType>
-void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m)
+void createRandomProjectionOfRank(int desired_rank, int rows, int cols, MatrixType& m)
 {
   typedef typename ei_traits<MatrixType>::Scalar Scalar;
   enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index 16eb27c52..abee32184 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -36,7 +36,7 @@ template<typename MatrixType> void qr()
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
   MatrixType m1;
-  createRandomMatrixOfRank(rank,rows,cols,m1);
+  createRandomProjectionOfRank(rank,rows,cols,m1);
   ColPivHouseholderQR<MatrixType> qr(m1);
   VERIFY_IS_APPROX(rank, qr.rank());
   VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
@@ -64,7 +64,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
   typedef typename MatrixType::Scalar Scalar;
   int rank = ei_random<int>(1, std::min(int(Rows), int(Cols))-1);
   Matrix<Scalar,Rows,Cols> m1;
-  createRandomMatrixOfRank(rank,Rows,Cols,m1);
+  createRandomProjectionOfRank(rank,Rows,Cols,m1);
   ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
   VERIFY_IS_APPROX(rank, qr.rank());
   VERIFY(Cols - qr.rank() == qr.dimensionOfKernel());
diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp
index c82ba4c7e..60255f94c 100644
--- a/test/qr_fullpivoting.cpp
+++ b/test/qr_fullpivoting.cpp
@@ -35,7 +35,7 @@ template<typename MatrixType> void qr()
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
   MatrixType m1;
-  createRandomMatrixOfRank(rank,rows,cols,m1);
+  createRandomProjectionOfRank(rank,rows,cols,m1);
   FullPivHouseholderQR<MatrixType> qr(m1);
   VERIFY_IS_APPROX(rank, qr.rank());
   VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
-- 
cgit v1.2.3


From d92df336ad21c7f8e0289f8ac3084b6313a17fe4 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Tue, 23 Feb 2010 15:40:24 -0500
Subject: Further LU test improvements. I'm not aware of any test failures
 anymore, not even with huge numbers of repetitions. Finally the
 createRandomMatrixOfRank() function is renamed to createRandomPIMatrixOfRank,
 where PI stands for 'partial isometry', that is, a matrix whose singular
 values are 0 or 1.

---
 Eigen/src/LU/FullPivLU.h |  7 ++++--
 test/inverse.cpp         |  2 +-
 test/lu.cpp              | 65 ++++++++++++++++++------------------------------
 test/main.h              | 12 ++++++---
 test/qr_colpivoting.cpp  |  4 +--
 test/qr_fullpivoting.cpp |  2 +-
 6 files changed, 42 insertions(+), 50 deletions(-)

(limited to 'Eigen/src/LU')

diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index ec551645b..0a305d52b 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -422,8 +422,11 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
     // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix
     if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();
 
-    // if the pivot (hence the corner) is exactly zero, terminate to avoid generating nan/inf values
-    if(ei_abs(biggest_in_corner) < cutoff)
+    // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values.
+    // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in
+    // their pseudo-code, results in numerical instability! The cutoff here has been validated
+    // by running the unit test 'lu' with many repetitions.
+    if(biggest_in_corner < cutoff)
     {
       // before exiting, make sure to initialize the still uninitialized transpositions
       // in a sane state without destroying what we already have.
diff --git a/test/inverse.cpp b/test/inverse.cpp
index 3f6138e0c..1e567ad14 100644
--- a/test/inverse.cpp
+++ b/test/inverse.cpp
@@ -42,7 +42,7 @@ template<typename MatrixType> void inverse(const MatrixType& m)
              m2(rows, cols),
              mzero = MatrixType::Zero(rows, cols),
              identity = MatrixType::Identity(rows, rows);
-  createRandomProjectionOfRank(rows,rows,rows,m1);
+  createRandomPIMatrixOfRank(rows,rows,rows,m1);
   m2 = m1.inverse();
   VERIFY_IS_APPROX(m1, m2.inverse() );
 
diff --git a/test/lu.cpp b/test/lu.cpp
index 02f6ec805..442202a33 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -28,9 +28,6 @@ using namespace std;
 
 template<typename MatrixType> void lu_non_invertible()
 {
-  static int times_called = 0;
-  times_called++;
-  
   typedef typename MatrixType::Scalar Scalar;
   typedef typename MatrixType::RealScalar RealScalar;
   /* this test covers the following files:
@@ -55,11 +52,16 @@ template<typename MatrixType> void lu_non_invertible()
     cols2 = cols = MatrixType::ColsAtCompileTime;
   }
 
+  enum {
+    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+    ColsAtCompileTime = MatrixType::ColsAtCompileTime
+  };
   typedef typename ei_kernel_retval_base<FullPivLU<MatrixType> >::ReturnType KernelMatrixType;
   typedef typename ei_image_retval_base<FullPivLU<MatrixType> >::ReturnType ImageMatrixType;
-  typedef Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> DynamicMatrixType;
-  typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime>
+  typedef Matrix<typename MatrixType::Scalar, ColsAtCompileTime, ColsAtCompileTime>
           CMatrixType;
+  typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, RowsAtCompileTime>
+          RMatrixType;
 
   int rank = ei_random<int>(1, std::min(rows, cols)-1);
 
@@ -68,26 +70,21 @@ template<typename MatrixType> void lu_non_invertible()
   
   MatrixType m1(rows, cols), m3(rows, cols2);
   CMatrixType m2(cols, cols2);
-  createRandomProjectionOfRank(rank, rows, cols, m1);
+  createRandomPIMatrixOfRank(rank, rows, cols, m1);
 
   FullPivLU<MatrixType> lu;
 
-  // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank of projections.
-  // So it's not clear at all the epsilon should play any role there.
+  // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank
+  // of singular values are either 0 or 1.
+  // So it's not clear at all that the epsilon should play any role there.
   lu.setThreshold(RealScalar(0.01));
   lu.compute(m1);
 
-  // FIXME need better way to construct trapezoid matrices. extend triangularView to support rectangular.
-  DynamicMatrixType u(rows,cols);
-  for(int i = 0; i < rows; i++)
-    for(int j = 0; j < cols; j++)
-      u(i,j) = i>j ? Scalar(0) : lu.matrixLU()(i,j);
-  DynamicMatrixType l(rows,rows);
-  for(int i = 0; i < rows; i++)
-    for(int j = 0; j < rows; j++)
-      l(i,j) = (i<j || j>=cols) ? Scalar(0)
-             : i==j ? Scalar(1)
-             : lu.matrixLU()(i,j);
+  MatrixType u(rows,cols);
+  u = lu.matrixLU().template triangularView<Upper>();
+  RMatrixType l = RMatrixType::Identity(rows,rows);
+  l.block(0,0,rows,std::min(rows,cols)).template triangularView<StrictlyLower>()
+    = lu.matrixLU().block(0,0,rows,std::min(rows,cols));
   
   VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u);
   
@@ -101,20 +98,8 @@ template<typename MatrixType> void lu_non_invertible()
   VERIFY(!lu.isSurjective());
   VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
   VERIFY(m1image.fullPivLu().rank() == rank);
+  VERIFY_IS_APPROX(m1 * m1.adjoint() * m1image, m1image);
 
-  // The following test is damn hard to get to succeed over a large number of repetitions.
-  // We're checking that the image is indeed the image, i.e. adding it as new columns doesn't increase the rank.
-  // Since we've already tested rank() above, the point here is not to test rank(), it is to test image().
-  // Since image() is implemented in a very simple way that doesn't leave much room for choice, the occasional
-  // errors that we get here (one in 1e+4 repetitions roughly) are probably just a sign that it's a really
-  // hard test, so we just limit how many times it's run.
-  if(times_called < 100)
-  {
-    DynamicMatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
-    sidebyside << m1, m1image;
-    VERIFY(sidebyside.fullPivLu().rank() == rank);
-  }
-  
   m2 = CMatrixType::Random(cols,cols2);
   m3 = m1*m2;
   m2 = CMatrixType::Random(cols,cols2);
@@ -128,20 +113,18 @@ template<typename MatrixType> void lu_invertible()
   /* this test covers the following files:
      LU.h
   */
+  typedef typename MatrixType::Scalar Scalar;
   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
   int size = ei_random<int>(1,200);
 
   MatrixType m1(size, size), m2(size, size), m3(size, size);
-  m1 = MatrixType::Random(size,size);
-
-  if (ei_is_same_type<RealScalar,float>::ret)
-  {
-    // let's build a matrix more stable to inverse
-    MatrixType a = MatrixType::Random(size,size*2);
-    m1 += a * a.adjoint();
-  }
+  FullPivLU<MatrixType> lu;
+  lu.setThreshold(0.01);
+  do {
+    m1 = MatrixType::Random(size,size);
+    lu.compute(m1);
+  } while(!lu.isInvertible());
 
-  FullPivLU<MatrixType> lu(m1);
   VERIFY(0 == lu.dimensionOfKernel());
   VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector
   VERIFY(size == lu.rank());
diff --git a/test/main.h b/test/main.h
index 6d296b2e3..96324de33 100644
--- a/test/main.h
+++ b/test/main.h
@@ -148,7 +148,7 @@ namespace Eigen
 
 #define EIGEN_INTERNAL_DEBUGGING
 #define EIGEN_NICE_RANDOM
-#include <Eigen/QR> // required for createRandomProjectionOfRank
+#include <Eigen/QR> // required for createRandomPIMatrixOfRank
 
 
 #define VERIFY(a) do { if (!(a)) { \
@@ -342,8 +342,13 @@ inline bool test_isUnitary(const MatrixBase<Derived>& m)
   return m.isUnitary(test_precision<typename ei_traits<Derived>::Scalar>());
 }
 
+/** Creates a random Partial Isometry matrix of given rank.
+  *
+  * A partial isometry is a matrix all of whose singular values are either 0 or 1.
+  * This is very useful to test rank-revealing algorithms.
+  */
 template<typename MatrixType>
-void createRandomProjectionOfRank(int desired_rank, int rows, int cols, MatrixType& m)
+void createRandomPIMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m)
 {
   typedef typename ei_traits<MatrixType>::Scalar Scalar;
   enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
@@ -360,7 +365,8 @@ void createRandomProjectionOfRank(int desired_rank, int rows, int cols, MatrixTy
 
   if(desired_rank == 1)
   {
-    m = VectorType::Random(rows) * VectorType::Random(cols).transpose();
+    // here we normalize the vectors to get a partial isometry
+    m = VectorType::Random(rows).normalized() * VectorType::Random(cols).normalized().transpose();
     return;
   }
 
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index abee32184..96cc66316 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -36,7 +36,7 @@ template<typename MatrixType> void qr()
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
   MatrixType m1;
-  createRandomProjectionOfRank(rank,rows,cols,m1);
+  createRandomPIMatrixOfRank(rank,rows,cols,m1);
   ColPivHouseholderQR<MatrixType> qr(m1);
   VERIFY_IS_APPROX(rank, qr.rank());
   VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
@@ -64,7 +64,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
   typedef typename MatrixType::Scalar Scalar;
   int rank = ei_random<int>(1, std::min(int(Rows), int(Cols))-1);
   Matrix<Scalar,Rows,Cols> m1;
-  createRandomProjectionOfRank(rank,Rows,Cols,m1);
+  createRandomPIMatrixOfRank(rank,Rows,Cols,m1);
   ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
   VERIFY_IS_APPROX(rank, qr.rank());
   VERIFY(Cols - qr.rank() == qr.dimensionOfKernel());
diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp
index 60255f94c..7ad3af1fe 100644
--- a/test/qr_fullpivoting.cpp
+++ b/test/qr_fullpivoting.cpp
@@ -35,7 +35,7 @@ template<typename MatrixType> void qr()
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
   MatrixType m1;
-  createRandomProjectionOfRank(rank,rows,cols,m1);
+  createRandomPIMatrixOfRank(rank,rows,cols,m1);
   FullPivHouseholderQR<MatrixType> qr(m1);
   VERIFY_IS_APPROX(rank, qr.rank());
   VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
-- 
cgit v1.2.3


From 7c98c04412322e56b3b6f7e235bc7ebb61ab6b43 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Wed, 24 Feb 2010 19:16:10 +0100
Subject: add reconstructedMatrix() to LLT, and LUs => they show that some
 improvements have still to be done    for permutations, tr*tr, trapezoidal
 matrices

---
 Eigen/src/Cholesky/LDLT.h   |  4 ++--
 Eigen/src/Cholesky/LLT.h    | 12 ++++++++++++
 Eigen/src/LU/FullPivLU.h    | 29 +++++++++++++++++++++++++++++
 Eigen/src/LU/PartialPivLU.h | 20 ++++++++++++++++++++
 test/cholesky.cpp           |  7 +++----
 test/lu.cpp                 | 21 +++++++++++++++++++++
 6 files changed, 87 insertions(+), 6 deletions(-)

(limited to 'Eigen/src/LU')

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 8cfc256bb..8699fe7e0 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -155,7 +155,7 @@ template<typename _MatrixType> class LDLT
       return m_matrix;
     }
 
-    const MatrixType reconstructedMatrix() const;
+    MatrixType reconstructedMatrix() const;
 
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
@@ -324,7 +324,7 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
  * i.e., it returns the product: P^T L D L^* P.
  * This function is provided for debug purpose. */
 template<typename MatrixType>
-const MatrixType LDLT<MatrixType>::reconstructedMatrix() const
+MatrixType LDLT<MatrixType>::reconstructedMatrix() const
 {
   ei_assert(m_isInitialized && "LDLT is not initialized.");
   const int size = m_matrix.rows();
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 96e1e5f73..2e8df7661 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -133,6 +133,8 @@ template<typename _MatrixType, int _UpLo> class LLT
       return m_matrix;
     }
 
+    MatrixType reconstructedMatrix() const;
+
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
 
@@ -295,6 +297,16 @@ bool LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
   return true;
 }
 
+/** \returns the matrix represented by the decomposition,
+ * i.e., it returns the product: L L^*.
+ * This function is provided for debug purpose. */
+template<typename MatrixType, int _UpLo>
+MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
+{
+  ei_assert(m_isInitialized && "LLT is not initialized.");
+  return matrixL() * matrixL().adjoint().toDenseMatrix();
+}
+
 /** \cholesky_module
   * \returns the LLT decomposition of \c *this
   */
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 0a305d52b..cd63b9ec7 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -361,6 +361,8 @@ template<typename _MatrixType> class FullPivLU
                (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
     }
 
+    MatrixType reconstructedMatrix() const;
+
     inline int rows() const { return m_lu.rows(); }
     inline int cols() const { return m_lu.cols(); }
 
@@ -487,6 +489,33 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons
   return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
 }
 
+/** \returns the matrix represented by the decomposition,
+ * i.e., it returns the product: P^{-1} L U Q^{-1}.
+ * This function is provided for debug purpose. */
+template<typename MatrixType>
+MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
+{
+  ei_assert(m_isInitialized && "LU is not initialized.");
+  const int smalldim = std::min(m_lu.rows(), m_lu.cols());
+  // LU
+  MatrixType res(m_lu.rows(),m_lu.cols());
+  // FIXME the .toDenseMatrix() should not be needed...
+  res = m_lu.corner(TopLeft,m_lu.rows(),smalldim)
+            .template triangularView<UnitLower>().toDenseMatrix()
+      * m_lu.corner(TopLeft,smalldim,m_lu.cols())
+            .template triangularView<Upper>().toDenseMatrix();
+
+  // P^{-1}(LU)
+  // FIXME implement inplace permutation
+  res = (m_p.inverse() * res).eval();
+
+  // (P^{-1}LU)Q^{-1}
+  // FIXME implement inplace permutation
+  res = (res * m_q.inverse()).eval();
+
+  return res;
+}
+
 /********* Implementation of kernel() **************************************************/
 
 template<typename _MatrixType>
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index ed2354d78..fcffc2458 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -165,6 +165,8 @@ template<typename _MatrixType> class PartialPivLU
       */
     typename ei_traits<MatrixType>::Scalar determinant() const;
 
+    MatrixType reconstructedMatrix() const;
+
     inline int rows() const { return m_lu.rows(); }
     inline int cols() const { return m_lu.cols(); }
 
@@ -400,6 +402,24 @@ typename ei_traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() c
   return Scalar(m_det_p) * m_lu.diagonal().prod();
 }
 
+/** \returns the matrix represented by the decomposition,
+ * i.e., it returns the product: P^{-1} L U.
+ * This function is provided for debug purpose. */
+template<typename MatrixType>
+MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
+{
+  ei_assert(m_isInitialized && "LU is not initialized.");
+  // LU
+  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
+                 * m_lu.template triangularView<Upper>();
+  
+  // P^{-1}(LU)
+  // FIXME implement inplace permutation
+  res = (m_p.inverse() * res).eval();
+
+  return res;
+}
+
 /***** Implementation of solve() *****************************************************/
 
 template<typename _MatrixType, typename Rhs>
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 1bb808d20..a446f5d73 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -95,7 +95,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
 
   {
     LLT<SquareMatrixType,Lower> chollo(symmLo);
-    VERIFY_IS_APPROX(symm, chollo.matrixL().toDenseMatrix() * chollo.matrixL().adjoint().toDenseMatrix());
+    VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
     vecX = chollo.solve(vecB);
     VERIFY_IS_APPROX(symm * vecX, vecB);
     matX = chollo.solve(matB);
@@ -103,7 +103,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
 
     // test the upper mode
     LLT<SquareMatrixType,Upper> cholup(symmUp);
-    VERIFY_IS_APPROX(symm, cholup.matrixL().toDenseMatrix() * cholup.matrixL().adjoint().toDenseMatrix());
+    VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix());
     vecX = cholup.solve(vecB);
     VERIFY_IS_APPROX(symm * vecX, vecB);
     matX = cholup.solve(matB);
@@ -119,8 +119,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
 
   {
     LDLT<SquareMatrixType> ldlt(symm);
-    // TODO(keir): This doesn't make sense now that LDLT pivots.
-    //VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint());
+    VERIFY_IS_APPROX(symm, ldlt.reconstructedMatrix());
     vecX = ldlt.solve(vecB);
     VERIFY_IS_APPROX(symm * vecX, vecB);
     matX = ldlt.solve(matB);
diff --git a/test/lu.cpp b/test/lu.cpp
index 442202a33..1ed38cb2b 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -91,6 +91,7 @@ template<typename MatrixType> void lu_non_invertible()
   KernelMatrixType m1kernel = lu.kernel();
   ImageMatrixType m1image = lu.image(m1);
 
+  VERIFY_IS_APPROX(m1, lu.reconstructedMatrix());
   VERIFY(rank == lu.rank());
   VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
   VERIFY(!lu.isInjective());
@@ -125,6 +126,7 @@ template<typename MatrixType> void lu_invertible()
     lu.compute(m1);
   } while(!lu.isInvertible());
 
+  VERIFY_IS_APPROX(m1, lu.reconstructedMatrix());
   VERIFY(0 == lu.dimensionOfKernel());
   VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector
   VERIFY(size == lu.rank());
@@ -138,6 +140,23 @@ template<typename MatrixType> void lu_invertible()
   VERIFY_IS_APPROX(m2, lu.inverse()*m3);
 }
 
+template<typename MatrixType> void lu_partial_piv()
+{
+  /* this test covers the following files:
+     PartialPivLU.h
+  */
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+  int rows = ei_random<int>(1,4);
+  int cols = rows;
+
+  MatrixType m1(cols, rows);
+  m1.setRandom();
+  PartialPivLU<MatrixType> plu(m1);
+
+  VERIFY_IS_APPROX(m1, plu.reconstructedMatrix());
+}
+
 template<typename MatrixType> void lu_verify_assert()
 {
   MatrixType tmp;
@@ -180,6 +199,7 @@ void test_lu()
     
     CALL_SUBTEST_4( lu_non_invertible<MatrixXd>() );
     CALL_SUBTEST_4( lu_invertible<MatrixXd>() );
+    CALL_SUBTEST_4( lu_partial_piv<MatrixXd>() );
     CALL_SUBTEST_4( lu_verify_assert<MatrixXd>() );
     
     CALL_SUBTEST_5( lu_non_invertible<MatrixXcf>() );
@@ -188,6 +208,7 @@ void test_lu()
     
     CALL_SUBTEST_6( lu_non_invertible<MatrixXcd>() );
     CALL_SUBTEST_6( lu_invertible<MatrixXcd>() );
+    CALL_SUBTEST_6( lu_partial_piv<MatrixXcd>() );
     CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() );
 
     CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() ));
-- 
cgit v1.2.3


From 959a1b5d6335833e9ad49a088502705bb6967ff5 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Thu, 25 Feb 2010 16:30:58 +0100
Subject: detect and implement inplace permutations

---
 Eigen/src/Core/PermutationMatrix.h | 53 ++++++++++++++++++++++++++++----------
 Eigen/src/Core/Transpose.h         | 19 --------------
 Eigen/src/Core/util/BlasUtil.h     | 18 +++++++++++++
 Eigen/src/LU/FullPivLU.h           |  8 +++---
 Eigen/src/LU/PartialPivLU.h        |  5 ++--
 test/permutationmatrices.cpp       | 19 +++++++++++++-
 6 files changed, 80 insertions(+), 42 deletions(-)

(limited to 'Eigen/src/LU')

diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h
index c42812ec8..46884dc3f 100644
--- a/Eigen/src/Core/PermutationMatrix.h
+++ b/Eigen/src/Core/PermutationMatrix.h
@@ -326,21 +326,46 @@ struct ei_permut_matrix_product_retval
     template<typename Dest> inline void evalTo(Dest& dst) const
     {
       const int n = Side==OnTheLeft ? rows() : cols();
-      for(int i = 0; i < n; ++i)
+
+      if(ei_is_same_type<MatrixTypeNestedCleaned,Dest>::ret && ei_extract_data(dst) == ei_extract_data(m_matrix))
+      {
+        // apply the permutation inplace
+        Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
+        mask.fill(false);
+        int r = 0;
+        while(r < m_permutation.size())
+        {
+          // search for the next seed
+          while(r<m_permutation.size() && mask[r]) r++;
+          if(r>=m_permutation.size())
+            break;
+          // we got one, let's follow it until we are back to the seed
+          int k0 = r++;
+          int kPrev = k0;
+          mask.coeffRef(k0) = true;
+          for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
+          {
+                  Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
+            .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
+                       (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
+
+            mask.coeffRef(k) = true;
+            kPrev = k;
+          }
+        }
+      }
+      else
       {
-        Block<
-          Dest,
-          Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime,
-          Side==OnTheRight ? 1 : Dest::ColsAtCompileTime
-        >(dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
-
-        =
-
-        Block<
-          MatrixTypeNestedCleaned,
-          Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,
-          Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime
-        >(m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
+        for(int i = 0; i < n; ++i)
+        {
+          Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
+               (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
+
+          =
+
+          Block<MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
+               (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
+        }
       }
     }
 
diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h
index bd06d8464..6c0e50de2 100644
--- a/Eigen/src/Core/Transpose.h
+++ b/Eigen/src/Core/Transpose.h
@@ -295,25 +295,6 @@ struct ei_blas_traits<SelfCwiseBinaryOp<BinOp,NestedXpr> >
   static inline const XprType extract(const XprType& x) { return x; }
 };
 
-
-template<typename T, int Access=ei_blas_traits<T>::ActualAccess>
-struct ei_extract_data_selector {
-  static typename T::Scalar* run(const T& m)
-  {
-    return &ei_blas_traits<T>::extract(m).const_cast_derived().coeffRef(0,0);
-  }
-};
-
-template<typename T>
-struct ei_extract_data_selector<T,NoDirectAccess> {
-  static typename T::Scalar* run(const T&) { return 0; }
-};
-
-template<typename T> typename T::Scalar* ei_extract_data(const T& m)
-{
-  return ei_extract_data_selector<T>::run(m);
-}
-
 template<typename Scalar, bool DestIsTranposed, typename OtherDerived>
 struct ei_check_transpose_aliasing_selector
 {
diff --git a/Eigen/src/Core/util/BlasUtil.h b/Eigen/src/Core/util/BlasUtil.h
index 2ca463d5d..4d216d77a 100644
--- a/Eigen/src/Core/util/BlasUtil.h
+++ b/Eigen/src/Core/util/BlasUtil.h
@@ -236,4 +236,22 @@ struct ei_blas_traits<Transpose<NestedXpr> >
   static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x.nestedExpression()); }
 };
 
+template<typename T, int Access=ei_blas_traits<T>::ActualAccess>
+struct ei_extract_data_selector {
+  static const typename T::Scalar* run(const T& m)
+  {
+    return &ei_blas_traits<T>::extract(m).const_cast_derived().coeffRef(0,0); // FIXME this should be .data()
+  }
+};
+
+template<typename T>
+struct ei_extract_data_selector<T,NoDirectAccess> {
+  static typename T::Scalar* run(const T&) { return 0; }
+};
+
+template<typename T> const typename T::Scalar* ei_extract_data(const T& m)
+{
+  return ei_extract_data_selector<T>::run(m);
+}
+
 #endif // EIGEN_BLASUTIL_H
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index cd63b9ec7..dea6ec41c 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -119,7 +119,7 @@ template<typename _MatrixType> class FullPivLU
       *          diagonal coefficient of U.
       */
     RealScalar maxPivot() const { return m_maxpivot; }
-    
+
     /** \returns the permutation matrix P
       *
       * \sa permutationQ()
@@ -506,12 +506,10 @@ MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
             .template triangularView<Upper>().toDenseMatrix();
 
   // P^{-1}(LU)
-  // FIXME implement inplace permutation
-  res = (m_p.inverse() * res).eval();
+  res = m_p.inverse() * res;
 
   // (P^{-1}LU)Q^{-1}
-  // FIXME implement inplace permutation
-  res = (res * m_q.inverse()).eval();
+  res = res * m_q.inverse();
 
   return res;
 }
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index fcffc2458..ad0d6b810 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -412,10 +412,9 @@ MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
   // LU
   MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
                  * m_lu.template triangularView<Upper>();
-  
+
   // P^{-1}(LU)
-  // FIXME implement inplace permutation
-  res = (m_p.inverse() * res).eval();
+  res = m_p.inverse() * res;
 
   return res;
 }
diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp
index ae1bd8b85..89142d910 100644
--- a/test/permutationmatrices.cpp
+++ b/test/permutationmatrices.cpp
@@ -86,6 +86,23 @@ template<typename MatrixType> void permutationmatrices(const MatrixType& m)
   identityp.setIdentity(rows);
   VERIFY_IS_APPROX(m_original, identityp*m_original);
 
+  // check inplace permutations
+  m_permuted = m_original;
+  m_permuted = lp.inverse() * m_permuted;
+  VERIFY_IS_APPROX(m_permuted, lp.inverse()*m_original);
+
+  m_permuted = m_original;
+  m_permuted = m_permuted * rp.inverse();
+  VERIFY_IS_APPROX(m_permuted, m_original*rp.inverse());
+
+  m_permuted = m_original;
+  m_permuted = lp * m_permuted;
+  VERIFY_IS_APPROX(m_permuted, lp*m_original);
+
+  m_permuted = m_original;
+  m_permuted = m_permuted * rp;
+  VERIFY_IS_APPROX(m_permuted, m_original*rp);
+
   if(rows>1 && cols>1)
   {
     lp2 = lp;
@@ -114,7 +131,7 @@ void test_permutationmatrices()
     CALL_SUBTEST_2( permutationmatrices(Matrix3f()) );
     CALL_SUBTEST_3( permutationmatrices(Matrix<double,3,3,RowMajor>()) );
     CALL_SUBTEST_4( permutationmatrices(Matrix4d()) );
-    CALL_SUBTEST_5( permutationmatrices(Matrix<double,4,6>()) );
+    CALL_SUBTEST_5( permutationmatrices(Matrix<double,40,60>()) );
     CALL_SUBTEST_6( permutationmatrices(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 30)) );
     CALL_SUBTEST_7( permutationmatrices(MatrixXcf(15, 10)) );
   }
-- 
cgit v1.2.3