From 1a1b2e9f27db619303e7f212f9bf5c58a2dd988c Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Fri, 10 Jul 2009 10:41:26 +0200
Subject: finally directly calling the low-level products is faster

---
 test/triangular.cpp | 54 ++++++++++++++++++++---------------------------------
 1 file changed, 20 insertions(+), 34 deletions(-)

(limited to 'test/triangular.cpp')

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)) );
   }
 }
-- 
cgit v1.2.3