src/LU/Determinant.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_DETERMINANT_H
0011 #define EIGEN_DETERMINANT_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 template<typename Derived>
0018 EIGEN_DEVICE_FUNC
0019 inline const typename Derived::Scalar bruteforce_det3_helper
0020 (const MatrixBase<Derived>& matrix, int a, int b, int c)
0021 {
0022   return matrix.coeff(0,a)
0023          * (matrix.coeff(1,b) * matrix.coeff(2,c) - matrix.coeff(1,c) * matrix.coeff(2,b));
0024 }
0025
0026 template<typename Derived,
0027          int DeterminantType = Derived::RowsAtCompileTime
0028 > struct determinant_impl
0029 {
0030   static inline typename traits<Derived>::Scalar run(const Derived& m)
0031   {
0032     if(Derived::ColsAtCompileTime==Dynamic && m.rows()==0)
0033       return typename traits<Derived>::Scalar(1);
0034     return m.partialPivLu().determinant();
0035   }
0036 };
0037
0038 template<typename Derived> struct determinant_impl<Derived, 1>
0039 {
0040   static inline EIGEN_DEVICE_FUNC
0041   typename traits<Derived>::Scalar run(const Derived& m)
0042   {
0043     return m.coeff(0,0);
0044   }
0045 };
0046
0047 template<typename Derived> struct determinant_impl<Derived, 2>
0048 {
0049   static inline EIGEN_DEVICE_FUNC
0050   typename traits<Derived>::Scalar run(const Derived& m)
0051   {
0052     return m.coeff(0,0) * m.coeff(1,1) - m.coeff(1,0) * m.coeff(0,1);
0053   }
0054 };
0055
0056 template<typename Derived> struct determinant_impl<Derived, 3>
0057 {
0058   static inline EIGEN_DEVICE_FUNC
0059   typename traits<Derived>::Scalar run(const Derived& m)
0060   {
0061     return bruteforce_det3_helper(m,0,1,2)
0062           - bruteforce_det3_helper(m,1,0,2)
0063           + bruteforce_det3_helper(m,2,0,1);
0064   }
0065 };
0066
0067 template<typename Derived> struct determinant_impl<Derived, 4>
0068 {
0069   typedef typename traits<Derived>::Scalar Scalar;
0070   static EIGEN_DEVICE_FUNC
0071   Scalar run(const Derived& m)
0072   {
0073     Scalar d2_01 = det2(m, 0, 1);
0074     Scalar d2_02 = det2(m, 0, 2);
0075     Scalar d2_03 = det2(m, 0, 3);
0076     Scalar d2_12 = det2(m, 1, 2);
0077     Scalar d2_13 = det2(m, 1, 3);
0078     Scalar d2_23 = det2(m, 2, 3);
0079     Scalar d3_0 = det3(m, 1,d2_23, 2,d2_13, 3,d2_12);
0080     Scalar d3_1 = det3(m, 0,d2_23, 2,d2_03, 3,d2_02);
0081     Scalar d3_2 = det3(m, 0,d2_13, 1,d2_03, 3,d2_01);
0082     Scalar d3_3 = det3(m, 0,d2_12, 1,d2_02, 2,d2_01);
0083     return internal::pmadd(-m(0,3),d3_0, m(1,3)*d3_1) +
0084            internal::pmadd(-m(2,3),d3_2, m(3,3)*d3_3);
0085   }
0086 protected:
0087   static EIGEN_DEVICE_FUNC
0088   Scalar det2(const Derived& m, Index i0, Index i1)
0089   {
0090     return m(i0,0) * m(i1,1) - m(i1,0) * m(i0,1);
0091   }
0092
0093   static EIGEN_DEVICE_FUNC
0094   Scalar det3(const Derived& m, Index i0, const Scalar& d0, Index i1, const Scalar& d1, Index i2, const Scalar& d2)
0095   {
0096     return internal::pmadd(m(i0,2), d0, internal::pmadd(-m(i1,2), d1, m(i2,2)*d2));
0097   }
0098 };
0099
0100 } // end namespace internal
0101
0102 /** \lu_module
0103   *
0104   * \returns the determinant of this matrix
0105   */
0106 template<typename Derived>
0107 EIGEN_DEVICE_FUNC
0108 inline typename internal::traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
0109 {
0110   eigen_assert(rows() == cols());
0111   typedef typename internal::nested_eval<Derived,Base::RowsAtCompileTime>::type Nested;
0112   return internal::determinant_impl<typename internal::remove_all<Nested>::type>::run(derived());
0113 }
0114
0115 } // end namespace Eigen
0116
0117 #endif // EIGEN_DETERMINANT_H