From 173e582e3ce7a7822cb40dbfe8855e642a5fc9f2 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Fri, 25 Apr 2008 23:13:20 +0000
Subject: [PATCH] added a tough test to check the determinant that currently
 fails

---
 test/CMakeLists.txt      |  1 +
 test/determinant.cpp     | 71 ++++++++++++++++++++++++++++++++++++++++
 test/linearstructure.cpp |  7 ++++
 test/main.h              |  1 +
 test/product.cpp         |  7 ++--
 5 files changed, 85 insertions(+), 2 deletions(-)
 create mode 100644 test/determinant.cpp

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 22bc2d23c..4a5ca235c 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -17,6 +17,7 @@ SET(test_SRCS
   smallvectors.cpp
   map.cpp
   cwiseop.cpp
+  determinant.cpp
 )
 QT4_AUTOMOC(${test_SRCS})
 
diff --git a/test/determinant.cpp b/test/determinant.cpp
new file mode 100644
index 000000000..974df58ef
--- /dev/null
+++ b/test/determinant.cpp
@@ -0,0 +1,71 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 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
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+#include <Eigen/LU>
+
+namespace Eigen {
+
+template<typename MatrixType> void nullDeterminant(const MatrixType& m)
+{
+  /* this test covers the following files:
+     Determinant.h
+  */
+  int rows = m.rows();
+  int cols = m.cols();
+
+  typedef typename MatrixType::Scalar Scalar;
+  typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> SquareMatrixType;
+  typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
+
+  MatrixType d(rows, cols);
+
+  // build a ill-conditionned matrix with a nul determinant
+  d.col(0).setOnes();
+  d.block(0,1, rows, cols-2).setRandom();
+  d.col(cols-1).setOnes();
+
+  for (int i=0 ; i<rows ; ++i)
+    d.row(i).block(0,1,1,cols-2) = d.row(i).block(0,1,1,cols-2).normalized();
+
+  SquareMatrixType covarianceMatrix = d.transpose() * d;
+
+//   std::cout << covarianceMatrix << "\n" << covarianceMatrix.determinant() << "\n";
+
+  VERIFY_IS_APPROX(covarianceMatrix.determinant(), Scalar(0));
+}
+
+void EigenTest::testDeterminant()
+{
+  for(int i = 0; i < m_repeat; i++) {
+    nullDeterminant(Matrix<float, 30, 3>());
+    nullDeterminant(Matrix<double, 30, 3>());
+    nullDeterminant(Matrix<float, 20, 4>());
+    nullDeterminant(Matrix<double, 20, 4>());
+//     nullDeterminant(MatrixXd(20,4));
+  }
+}
+
+} // namespace Eigen
diff --git a/test/linearstructure.cpp b/test/linearstructure.cpp
index 6e1364ec4..25d0bdc2c 100644
--- a/test/linearstructure.cpp
+++ b/test/linearstructure.cpp
@@ -89,6 +89,12 @@ template<typename MatrixType> void linearStructure(const MatrixType& m)
   VERIFY_IS_APPROX((m1*s1)(r,c), (m1(r,c))*s1);
   if(NumTraits<Scalar>::HasFloatingPoint)
     VERIFY_IS_APPROX((m1/s1)(r,c), (m1(r,c))/s1);
+
+  // use .block to disable vectorization and compare to the vectorized version
+  VERIFY_IS_APPROX(m1+m1.block(0,0,rows,cols), m1+m1);
+  VERIFY_IS_APPROX(m1.cwiseProduct(m1.block(0,0,rows,cols)), m1.cwiseProduct(m1));
+  VERIFY_IS_APPROX(m1 - m1.block(0,0,rows,cols), m1 - m1);
+  VERIFY_IS_APPROX(m1.block(0,0,rows,cols) * s1, m1 * s1);
 }
 
 void EigenTest::testLinearStructure()
@@ -97,6 +103,7 @@ void EigenTest::testLinearStructure()
     linearStructure(Matrix<float, 1, 1>());
     linearStructure(Matrix4d());
     linearStructure(MatrixXcf(3, 3));
+    linearStructure(MatrixXf(8, 12));
     linearStructure(MatrixXi(8, 12));
     linearStructure(MatrixXcd(20, 20));
   }
diff --git a/test/main.h b/test/main.h
index 67b5b456f..3aae7bc65 100644
--- a/test/main.h
+++ b/test/main.h
@@ -213,6 +213,7 @@ class EigenTest : public QObject
     void testSmallVectors();
     void testMap();
     void testCwiseops();
+    void testDeterminant();
   protected:
     int m_repeat;
 };
diff --git a/test/product.cpp b/test/product.cpp
index a6b6a537d..782be8f82 100644
--- a/test/product.cpp
+++ b/test/product.cpp
@@ -34,6 +34,7 @@ template<typename MatrixType> void product(const MatrixType& m)
 
   typedef typename MatrixType::Scalar Scalar;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
 
   int rows = m.rows();
   int cols = m.cols();
@@ -43,7 +44,8 @@ template<typename MatrixType> void product(const MatrixType& m)
   MatrixType m1 = MatrixType::random(rows, cols),
              m2 = MatrixType::random(rows, cols),
              m3(rows, cols),
-             mzero = MatrixType::zero(rows, cols),
+             mzero = MatrixType::zero(rows, cols);
+  SquareMatrixType
              identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
                               ::identity(rows, rows),
              square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
@@ -95,7 +97,8 @@ void EigenTest::testProduct()
     product(Matrix<float, 1, 1>());
     product(Matrix4d());
     product(MatrixXcf(3, 3));
-    product(MatrixXi(8, 12));
+    product(MatrixXf(13, 25));
+    product(MatrixXi(4, 4));
     product(MatrixXcd(20, 20));
   }