- 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

1 files changed, 68 insertions, 11 deletions

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