finally directly calling the low-level products is faster

4 files changed, 40 insertions, 90 deletions

diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 253f67282..4b3952a62 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -142,18 +142,18 @@ template<typename MatrixType> void cholesky_verify_assert()
 void test_cholesky()
 {
   for(int i = 0; i < g_repeat; i++) {
-//     CALL_SUBTEST( cholesky(Matrix<double,1,1>()) );
-//     CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
-//     CALL_SUBTEST( cholesky(Matrix2d()) );
-//     CALL_SUBTEST( cholesky(Matrix3f()) );
+    CALL_SUBTEST( cholesky(Matrix<double,1,1>()) );
+    CALL_SUBTEST( cholesky(MatrixXd(1,1)) );
+    CALL_SUBTEST( cholesky(Matrix2d()) );
+    CALL_SUBTEST( cholesky(Matrix3f()) );
     CALL_SUBTEST( cholesky(Matrix4d()) );
     CALL_SUBTEST( cholesky(MatrixXcd(7,7)) );
-//     CALL_SUBTEST( cholesky(MatrixXd(17,17)) );
-//     CALL_SUBTEST( cholesky(MatrixXf(200,200)) );
+    CALL_SUBTEST( cholesky(MatrixXd(17,17)) );
+    CALL_SUBTEST( cholesky(MatrixXf(200,200)) );
   }
 
-//   CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
-//   CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
-//   CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
-//   CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
+  CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
+  CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
+  CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
+  CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
 }
diff --git a/test/main.h b/test/main.h
index a6a780603..76572d9b3 100644
--- a/test/main.h
+++ b/test/main.h
@@ -28,6 +28,7 @@
 #include <iostream>
 #include <string>
 #include <vector>
+#include <typeinfo>
 
 #ifndef EIGEN_TEST_FUNC
 #error EIGEN_TEST_FUNC must be defined
diff --git a/test/product_extra.cpp b/test/product_extra.cpp
index d73974886..e750be65e 100644
--- a/test/product_extra.cpp
+++ b/test/product_extra.cpp
@@ -47,8 +47,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
              square2 = MatrixType::Random(cols, cols),
              res2 = MatrixType::Random(cols, cols);
   RowVectorType v1 = RowVectorType::Random(rows), vrres(rows);
-//              v2 = RowVectorType::Random(rows),
-//              vzero = RowVectorType::Zero(rows);
   ColVectorType vc2 = ColVectorType::Random(cols), vcres(cols);
   OtherMajorMatrixType tm1 = m1;
 
@@ -56,14 +54,14 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
          s2 = ei_random<Scalar>(),
          s3 = ei_random<Scalar>();
 
-  int c0 = ei_random<int>(0,cols/2-1),
-      c1 = ei_random<int>(cols/2,cols),
-      r0 = ei_random<int>(0,rows/2-1),
-      r1 = ei_random<int>(rows/2,rows);
+//   int c0 = ei_random<int>(0,cols/2-1),
+//       c1 = ei_random<int>(cols/2,cols),
+//       r0 = ei_random<int>(0,rows/2-1),
+//       r1 = ei_random<int>(rows/2,rows);
 
   // all the expressions in this test should be compiled as a single matrix product
   // TODO: add internal checks to verify that
-/*
+
   VERIFY_IS_APPROX(m1 * m2.adjoint(),  m1 * m2.adjoint().eval());
   VERIFY_IS_APPROX(m1.adjoint() * square.adjoint(),  m1.adjoint().eval() * square.adjoint().eval());
   VERIFY_IS_APPROX(m1.adjoint() * m2,  m1.adjoint().eval() * m2);
@@ -111,49 +109,14 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
 
   VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()),
                    (-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval());
-  */
-  // test with sub matrices
-  m2 = m1;
-  m3 = m1;
-
-//   std::cerr << (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).rows() << " " << (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).cols() << " == " << vrres.segment(r0,r1-r0).rows() << " " << vrres.segment(r0,r1-r0).cols() << "\n";
-//   m2.col(c0).segment(0,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).lazy();
-//   m3.col(c0).segment(0,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).eval();
-  Matrix<Scalar,Dynamic,1> a = m2.col(c0), b = a;
-  a.segment(5,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).lazy();
-  b.segment(5,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).eval();
-
-//   m2.col(c0).segment(0,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).lazy();
-//   m3.col(c0).segment(0,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).eval();
-//   if (!m2.isApprox(m3))
-  std::cerr << (a-b).cwise().abs().maxCoeff() << "\n";
-  VERIFY_IS_APPROX(a,b);
-//   VERIFY_IS_APPROX( vrres.segment(0,r1-r0).transpose().eval(),
-//                     v1.segment(0,r1-r0).transpose() + m1.block(r0,c0, r1-r0, c1-c0).eval() * (vc2.segment(c0,c1-c0)).eval());
+
 }
 
 void test_product_extra()
 {
   for(int i = 0; i < g_repeat; i++) {
-    int rows = ei_random<int>(2,10);
-    int cols = ei_random<int>(2,10);
-    int c0 = ei_random<int>(0,cols/2-1),
-        c1 = ei_random<int>(cols/2,cols),
-        r0 = ei_random<int>(0,rows/2-1),
-        r1 = ei_random<int>(rows/2,rows);
-
-    MatrixXf m1 = MatrixXf::Random(rows,cols), m2 = m1;
-    Matrix<float,Dynamic,1> a = m2.col(c0), b = a;
-    Matrix<float,Dynamic,1> vc2 = Matrix<float,Dynamic,1>::Random(cols);
-    if (1+r1-r0<rows) {
-      a.segment(1,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).lazy();
-      b.segment(1,r1-r0) += (m1.block(r0,c0, r1-r0, c1-c0) * vc2.segment(c0,c1-c0)).eval();
-      VERIFY_IS_APPROX(a,b);
-    }
-//     CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) );
-//     CALL_SUBTEST( product_extra(MatrixXd(ei_random<int>(1,320), ei_random<int>(1,320))) );
-//     CALL_SUBTEST( product(MatrixXi(ei_random<int>(1,320), ei_random<int>(1,320))) );
-//     CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) );
-//     CALL_SUBTEST( product(Matrix<float,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,320), ei_random<int>(1,320))) );
+    CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) );
+    CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) );
+    CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,50), ei_random<int>(1,50))) );
   }
 }
diff --git a/test/triangular.cpp b/test/triangular.cpp
index 3550d1a74..0c03e987e 100644
--- a/test/triangular.cpp
+++ b/test/triangular.cpp
@@ -1,7 +1,7 @@
 // This file is triangularView of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@gmail.com>
+// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -81,44 +81,35 @@ template<typename MatrixType> void triangular(const MatrixType& m)
   m1.template triangularView<Eigen::LowerTriangular>() = (m2.transpose() * m2).lazy();
   VERIFY_IS_APPROX(m3.template triangularView<Eigen::LowerTriangular>().toDense(), m1);
 
-  // VERIFY_IS_APPROX(m3.template triangularView<DiagonalBits>(), m3.diagonal().asDiagonal());
-
   m1 = MatrixType::Random(rows, cols);
   for (int i=0; i<rows; ++i)
     while (ei_abs2(m1(i,i))<1e-3) m1(i,i) = ei_random<Scalar>();
 
   Transpose<MatrixType> trm4(m4);
   // test back and forward subsitution
+  m3 = m1.template triangularView<Eigen::UpperTriangular>();
+  VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Eigen::LowerTriangular>().solve(m2)), largerEps));
+  m3 = m1.template triangularView<Eigen::LowerTriangular>();
+  VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Eigen::UpperTriangular>().solve(m2)), largerEps));
+  m3 = m1.template triangularView<Eigen::UpperTriangular>();
+  VERIFY(m2.isApprox(m3 * (m1.template triangularView<Eigen::UpperTriangular>().solve(m2)), largerEps));
   m3 = m1.template triangularView<Eigen::LowerTriangular>();
-//   VERIFY(m3.template triangularView<Eigen::LowerTriangular>().solve(m3).cwise().abs().isIdentity(test_precision<RealScalar>()));
-//   VERIFY(m3.transpose().template triangularView<Eigen::UpperTriangular>()
-//     .solve(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>()));
+  VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Eigen::LowerTriangular>().solve(m2)), largerEps));
+
   // check M * inv(L) using in place API
   m4 = m3;
   m3.transpose().template triangularView<Eigen::UpperTriangular>().solveInPlace(trm4);
-//   VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>()));
+  VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>()));
 
-//   m3 = m1.template triangularView<Eigen::UpperTriangular>();
-//   VERIFY(m3.template triangularView<Eigen::UpperTriangular>().solve(m3).cwise().abs().isIdentity(test_precision<RealScalar>()));
-//   VERIFY(m3.transpose().template triangularView<Eigen::LowerTriangular>()
-//     .solve(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>()));
-//   // check M * inv(U) using in place API
+  // check M * inv(U) using in place API
+  m3 = m1.template triangularView<Eigen::UpperTriangular>();
   m4 = m3;
   m3.transpose().template triangularView<Eigen::LowerTriangular>().solveInPlace(trm4);
   VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>()));
 
-//   m3 = m1.template triangularView<Eigen::UpperTriangular>();
-//   VERIFY(m2.isApprox(m3 * (m3.template triangularView<Eigen::UpperTriangular>().solve(m2)), largerEps));
-//   m3 = m1.template triangularView<Eigen::LowerTriangular>();
-
-//   std::cerr << (m2 -
-//             (m3 * (m3.template triangularView<Eigen::LowerTriangular>().solve(m2)))).cwise().abs() /*.maxCoeff()*/ << "\n\n";
-
-//   VERIFY(m2.isApprox(m3 * (m3.template triangularView<Eigen::LowerTriangular>().solve(m2)), largerEps));
-
   // check solve with unit diagonal
-//   m3 = m1.template triangularView<Eigen::UnitUpperTriangular>();
-//   VERIFY(m2.isApprox(m3 * (m1.template triangularView<Eigen::UnitUpperTriangular>().solve(m2)), largerEps));
+  m3 = m1.template triangularView<Eigen::UnitUpperTriangular>();
+  VERIFY(m2.isApprox(m3 * (m1.template triangularView<Eigen::UnitUpperTriangular>().solve(m2)), largerEps));
 
 //   VERIFY((  m1.template triangularView<Eigen::UpperTriangular>()
 //           * m2.template triangularView<Eigen::UpperTriangular>()).isUpperTriangular());
@@ -136,17 +127,12 @@ template<typename MatrixType> void triangular(const MatrixType& m)
 void test_triangular()
 {
   for(int i = 0; i < g_repeat ; i++) {
-//     CALL_SUBTEST( triangular(Matrix<float, 1, 1>()) );
-//     CALL_SUBTEST( triangular(Matrix<float, 2, 2>()) );
-//     CALL_SUBTEST( triangular(Matrix3d()) );
-//     CALL_SUBTEST( triangular(MatrixXcf(4, 4)) );
-//     CALL_SUBTEST( triangular(Matrix<std::complex<float>,8, 8>()) );
-//     CALL_SUBTEST( triangular(MatrixXd(1,1)) );
-//     CALL_SUBTEST( triangular(MatrixXd(2,2)) );
-//     CALL_SUBTEST( triangular(MatrixXd(3,3)) );
-//     CALL_SUBTEST( triangular(MatrixXd(5,5)) );
-//     CALL_SUBTEST( triangular(MatrixXd(8,8)) );
+    CALL_SUBTEST( triangular(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST( triangular(Matrix<float, 2, 2>()) );
+    CALL_SUBTEST( triangular(Matrix3d()) );
+    CALL_SUBTEST( triangular(MatrixXcf(4, 4)) );
+    CALL_SUBTEST( triangular(Matrix<std::complex<float>,8, 8>()) );
     CALL_SUBTEST( triangular(MatrixXd(17,17)) );
-//     CALL_SUBTEST( triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(5, 5)) );
+    CALL_SUBTEST( triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(5, 5)) );
   }
 }