From 5f350448e163c42b400c1191d790e3d2fa12d235 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Thu, 7 Aug 2008 21:48:21 +0000
Subject: - add kernel computation using the triangular solver - take advantage
 of the fact that our LU dec sorts the eigenvalues of U in decreasing order -
 add meta selector in determinant

---
 Eigen/src/LU/Determinant.h | 79 +++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 68 insertions(+), 11 deletions(-)

(limited to 'Eigen/src/LU/Determinant.h')

diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h
index 96e400564..f66836283 100644
--- a/Eigen/src/LU/Determinant.h
+++ b/Eigen/src/LU/Determinant.h
@@ -67,25 +67,82 @@ const typename Derived::Scalar ei_bruteforce_det(const MatrixBase<Derived>& m)
   }
 }
 
-/** \lu_module
-  *
-  * \returns the determinant of this matrix
-  */
-template<typename Derived>
-typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
+const int TriangularDeterminant = 0;
+
+template<typename Derived,
+         int DeterminantType =
+           (Derived::Flags & (UpperTriangularBit | LowerTriangularBit))
+           ? TriangularDeterminant : Derived::RowsAtCompileTime
+> struct ei_determinant_impl
 {
-  assert(rows() == cols());
-  if (Derived::Flags & (UpperTriangularBit | LowerTriangularBit))
+  static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
+  {
+    return m.lu().determinant();
+  }
+};
+
+template<typename Derived> struct ei_determinant_impl<Derived, TriangularDeterminant>
+{
+  static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
   {
     if (Derived::Flags & UnitDiagBit)
       return 1;
     else if (Derived::Flags & ZeroDiagBit)
       return 0;
     else
-      return derived().diagonal().redux(ei_scalar_product_op<Scalar>());
+      return m.diagonal().redux(ei_scalar_product_op<typename ei_traits<Derived>::Scalar>());
   }
-  else if(rows() <= 4) return ei_bruteforce_det(derived());
-  else return lu().determinant();
+};
+
+template<typename Derived> struct ei_determinant_impl<Derived, 1>
+{
+  static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
+  {
+    return m.coeff(0,0);
+  }
+};
+
+template<typename Derived> struct ei_determinant_impl<Derived, 2>
+{
+  static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
+  {
+    return m.coeff(0,0) * m.coeff(1,1) - m.coeff(1,0) * m.coeff(0,1);
+  }
+};
+
+template<typename Derived> struct ei_determinant_impl<Derived, 3>
+{
+  static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
+  {
+    return ei_bruteforce_det3_helper(m,0,1,2)
+          - ei_bruteforce_det3_helper(m,1,0,2)
+          + ei_bruteforce_det3_helper(m,2,0,1);
+  }
+};
+
+template<typename Derived> struct ei_determinant_impl<Derived, 4>
+{
+  static inline typename ei_traits<Derived>::Scalar run(const Derived& m)
+  {
+    // trick by Martin Costabel to compute 4x4 det with only 30 muls
+    return ei_bruteforce_det4_helper(m,0,1,2,3)
+          - ei_bruteforce_det4_helper(m,0,2,1,3)
+          + ei_bruteforce_det4_helper(m,0,3,1,2)
+          + ei_bruteforce_det4_helper(m,1,2,0,3)
+          - ei_bruteforce_det4_helper(m,1,3,0,2)
+          + ei_bruteforce_det4_helper(m,2,3,0,1);
+  }
+};
+
+/** \lu_module
+  *
+  * \returns the determinant of this matrix
+  */
+template<typename Derived>
+typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
+{
+  assert(rows() == cols());
+  return ei_determinant_impl<Derived>::run(derived());
 }
 
 #endif // EIGEN_DETERMINANT_H
-- 
cgit v1.2.3