Determinant.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_DETERMINANT_H
00011 #define EIGEN_DETERMINANT_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 template<typename Derived>
00018 inline const typename Derived::Scalar bruteforce_det3_helper
00019 (const MatrixBase<Derived>& matrix, int a, int b, int c)
00020 {
00021   return matrix.coeff(0,a)
00022          * (matrix.coeff(1,b) * matrix.coeff(2,c) - matrix.coeff(1,c) * matrix.coeff(2,b));
00023 }
00024 
00025 template<typename Derived>
00026 const typename Derived::Scalar bruteforce_det4_helper
00027 (const MatrixBase<Derived>& matrix, int j, int k, int m, int n)
00028 {
00029   return (matrix.coeff(j,0) * matrix.coeff(k,1) - matrix.coeff(k,0) * matrix.coeff(j,1))
00030        * (matrix.coeff(m,2) * matrix.coeff(n,3) - matrix.coeff(n,2) * matrix.coeff(m,3));
00031 }
00032 
00033 template<typename Derived,
00034          int DeterminantType = Derived::RowsAtCompileTime
00035 > struct determinant_impl
00036 {
00037   static inline typename traits<Derived>::Scalar run(const Derived& m)
00038   {
00039     if(Derived::ColsAtCompileTime==Dynamic && m.rows()==0)
00040       return typename traits<Derived>::Scalar(1);
00041     return m.partialPivLu().determinant();
00042   }
00043 };
00044 
00045 template<typename Derived> struct determinant_impl<Derived, 1>
00046 {
00047   static inline typename traits<Derived>::Scalar run(const Derived& m)
00048   {
00049     return m.coeff(0,0);
00050   }
00051 };
00052 
00053 template<typename Derived> struct determinant_impl<Derived, 2>
00054 {
00055   static inline typename traits<Derived>::Scalar run(const Derived& m)
00056   {
00057     return m.coeff(0,0) * m.coeff(1,1) - m.coeff(1,0) * m.coeff(0,1);
00058   }
00059 };
00060 
00061 template<typename Derived> struct determinant_impl<Derived, 3>
00062 {
00063   static inline typename traits<Derived>::Scalar run(const Derived& m)
00064   {
00065     return bruteforce_det3_helper(m,0,1,2)
00066           - bruteforce_det3_helper(m,1,0,2)
00067           + bruteforce_det3_helper(m,2,0,1);
00068   }
00069 };
00070 
00071 template<typename Derived> struct determinant_impl<Derived, 4>
00072 {
00073   static typename traits<Derived>::Scalar run(const Derived& m)
00074   {
00075     // trick by Martin Costabel to compute 4x4 det with only 30 muls
00076     return bruteforce_det4_helper(m,0,1,2,3)
00077           - bruteforce_det4_helper(m,0,2,1,3)
00078           + bruteforce_det4_helper(m,0,3,1,2)
00079           + bruteforce_det4_helper(m,1,2,0,3)
00080           - bruteforce_det4_helper(m,1,3,0,2)
00081           + bruteforce_det4_helper(m,2,3,0,1);
00082   }
00083 };
00084 
00085 } // end namespace internal
00086 
00087 /** \lu_module
00088   *
00089   * \returns the determinant of this matrix
00090   */
00091 template<typename Derived>
00092 inline typename internal::traits<Derived>::Scalar MatrixBase<Derived>::determinant () const
00093 {
00094   eigen_assert(rows() == cols());
00095   typedef typename internal::nested<Derived,Base::RowsAtCompileTime>::type Nested;
00096   return internal::determinant_impl<typename internal::remove_all<Nested>::type>::run(derived());
00097 }
00098 
00099 } // end namespace Eigen
00100 
00101 #endif // EIGEN_DETERMINANT_H