Inverse.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2010 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_INVERSE_H
00011 #define EIGEN_INVERSE_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 /**********************************
00018 *** General case implementation ***
00019 **********************************/
00020 
00021 template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
00022 struct compute_inverse
00023 {
00024   static inline void run(const MatrixType& matrix, ResultType& result)
00025   {
00026     result = matrix.partialPivLu().inverse();
00027   }
00028 };
00029 
00030 template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
00031 struct compute_inverse_and_det_with_check { /* nothing! general case not supported. */ };
00032 
00033 /****************************
00034 *** Size 1 implementation ***
00035 ****************************/
00036 
00037 template<typename MatrixType, typename ResultType>
00038 struct compute_inverse<MatrixType, ResultType, 1>
00039 {
00040   static inline void run(const MatrixType& matrix, ResultType& result)
00041   {
00042     typedef typename MatrixType::Scalar Scalar;
00043     result.coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
00044   }
00045 };
00046 
00047 template<typename MatrixType, typename ResultType>
00048 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 1>
00049 {
00050   static inline void run(
00051     const MatrixType& matrix,
00052     const typename MatrixType::RealScalar& absDeterminantThreshold,
00053     ResultType& result,
00054     typename ResultType::Scalar& determinant,
00055     bool& invertible
00056   )
00057   {
00058     determinant = matrix.coeff(0,0);
00059     invertible = abs(determinant) > absDeterminantThreshold;
00060     if(invertible) result.coeffRef(0,0) = typename ResultType::Scalar(1) / determinant;
00061   }
00062 };
00063 
00064 /****************************
00065 *** Size 2 implementation ***
00066 ****************************/
00067 
00068 template<typename MatrixType, typename ResultType>
00069 inline void compute_inverse_size2_helper(
00070     const MatrixType& matrix, const typename ResultType::Scalar& invdet,
00071     ResultType& result)
00072 {
00073   result.coeffRef(0,0) = matrix.coeff(1,1) * invdet;
00074   result.coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
00075   result.coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
00076   result.coeffRef(1,1) = matrix.coeff(0,0) * invdet;
00077 }
00078 
00079 template<typename MatrixType, typename ResultType>
00080 struct compute_inverse<MatrixType, ResultType, 2>
00081 {
00082   static inline void run(const MatrixType& matrix, ResultType& result)
00083   {
00084     typedef typename ResultType::Scalar Scalar;
00085     const Scalar invdet = typename MatrixType::Scalar(1) / matrix.determinant();
00086     compute_inverse_size2_helper(matrix, invdet, result);
00087   }
00088 };
00089 
00090 template<typename MatrixType, typename ResultType>
00091 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 2>
00092 {
00093   static inline void run(
00094     const MatrixType& matrix,
00095     const typename MatrixType::RealScalar& absDeterminantThreshold,
00096     ResultType& inverse,
00097     typename ResultType::Scalar& determinant,
00098     bool& invertible
00099   )
00100   {
00101     typedef typename ResultType::Scalar Scalar;
00102     determinant = matrix.determinant();
00103     invertible = abs(determinant) > absDeterminantThreshold;
00104     if(!invertible) return;
00105     const Scalar invdet = Scalar(1) / determinant;
00106     compute_inverse_size2_helper(matrix, invdet, inverse);
00107   }
00108 };
00109 
00110 /****************************
00111 *** Size 3 implementation ***
00112 ****************************/
00113 
00114 template<typename MatrixType, int i, int j>
00115 inline typename MatrixType::Scalar cofactor_3x3(const MatrixType& m)
00116 {
00117   enum {
00118     i1 = (i+1) % 3,
00119     i2 = (i+2) % 3,
00120     j1 = (j+1) % 3,
00121     j2 = (j+2) % 3
00122   };
00123   return m.coeff(i1, j1) * m.coeff(i2, j2)
00124        - m.coeff(i1, j2) * m.coeff(i2, j1);
00125 }
00126 
00127 template<typename MatrixType, typename ResultType>
00128 inline void compute_inverse_size3_helper(
00129     const MatrixType& matrix,
00130     const typename ResultType::Scalar& invdet,
00131     const Matrix<typename ResultType::Scalar,3,1>& cofactors_col0,
00132     ResultType& result)
00133 {
00134   result.row(0) = cofactors_col0 * invdet;
00135   result.coeffRef(1,0) =  cofactor_3x3<MatrixType,0,1>(matrix) * invdet;
00136   result.coeffRef(1,1) =  cofactor_3x3<MatrixType,1,1>(matrix) * invdet;
00137   result.coeffRef(1,2) =  cofactor_3x3<MatrixType,2,1>(matrix) * invdet;
00138   result.coeffRef(2,0) =  cofactor_3x3<MatrixType,0,2>(matrix) * invdet;
00139   result.coeffRef(2,1) =  cofactor_3x3<MatrixType,1,2>(matrix) * invdet;
00140   result.coeffRef(2,2) =  cofactor_3x3<MatrixType,2,2>(matrix) * invdet;
00141 }
00142 
00143 template<typename MatrixType, typename ResultType>
00144 struct compute_inverse<MatrixType, ResultType, 3>
00145 {
00146   static inline void run(const MatrixType& matrix, ResultType& result)
00147   {
00148     typedef typename ResultType::Scalar Scalar;
00149     Matrix<typename MatrixType::Scalar,3,1> cofactors_col0;
00150     cofactors_col0.coeffRef(0) =  cofactor_3x3<MatrixType,0,0>(matrix);
00151     cofactors_col0.coeffRef(1) =  cofactor_3x3<MatrixType,1,0>(matrix);
00152     cofactors_col0.coeffRef(2) =  cofactor_3x3<MatrixType,2,0>(matrix);
00153     const Scalar det = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
00154     const Scalar invdet = Scalar(1) / det;
00155     compute_inverse_size3_helper(matrix, invdet, cofactors_col0, result);
00156   }
00157 };
00158 
00159 template<typename MatrixType, typename ResultType>
00160 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 3>
00161 {
00162   static inline void run(
00163     const MatrixType& matrix,
00164     const typename MatrixType::RealScalar& absDeterminantThreshold,
00165     ResultType& inverse,
00166     typename ResultType::Scalar& determinant,
00167     bool& invertible
00168   )
00169   {
00170     typedef typename ResultType::Scalar Scalar;
00171     Matrix<Scalar,3,1> cofactors_col0;
00172     cofactors_col0.coeffRef(0) =  cofactor_3x3<MatrixType,0,0>(matrix);
00173     cofactors_col0.coeffRef(1) =  cofactor_3x3<MatrixType,1,0>(matrix);
00174     cofactors_col0.coeffRef(2) =  cofactor_3x3<MatrixType,2,0>(matrix);
00175     determinant = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
00176     invertible = abs(determinant) > absDeterminantThreshold;
00177     if(!invertible) return;
00178     const Scalar invdet = Scalar(1) / determinant;
00179     compute_inverse_size3_helper(matrix, invdet, cofactors_col0, inverse);
00180   }
00181 };
00182 
00183 /****************************
00184 *** Size 4 implementation ***
00185 ****************************/
00186 
00187 template<typename Derived>
00188 inline const typename Derived::Scalar general_det3_helper
00189 (const MatrixBase<Derived>& matrix, int i1, int i2, int i3, int j1, int j2, int j3)
00190 {
00191   return matrix.coeff(i1,j1)
00192          * (matrix.coeff(i2,j2) * matrix.coeff(i3,j3) - matrix.coeff(i2,j3) * matrix.coeff(i3,j2));
00193 }
00194 
00195 template<typename MatrixType, int i, int j>
00196 inline typename MatrixType::Scalar cofactor_4x4(const MatrixType& matrix)
00197 {
00198   enum {
00199     i1 = (i+1) % 4,
00200     i2 = (i+2) % 4,
00201     i3 = (i+3) % 4,
00202     j1 = (j+1) % 4,
00203     j2 = (j+2) % 4,
00204     j3 = (j+3) % 4
00205   };
00206   return general_det3_helper(matrix, i1, i2, i3, j1, j2, j3)
00207        + general_det3_helper(matrix, i2, i3, i1, j1, j2, j3)
00208        + general_det3_helper(matrix, i3, i1, i2, j1, j2, j3);
00209 }
00210 
00211 template<int Arch, typename Scalar, typename MatrixType, typename ResultType>
00212 struct compute_inverse_size4
00213 {
00214   static void run(const MatrixType& matrix, ResultType& result)
00215   {
00216     result.coeffRef(0,0) =  cofactor_4x4<MatrixType,0,0>(matrix);
00217     result.coeffRef(1,0) = -cofactor_4x4<MatrixType,0,1>(matrix);
00218     result.coeffRef(2,0) =  cofactor_4x4<MatrixType,0,2>(matrix);
00219     result.coeffRef(3,0) = -cofactor_4x4<MatrixType,0,3>(matrix);
00220     result.coeffRef(0,2) =  cofactor_4x4<MatrixType,2,0>(matrix);
00221     result.coeffRef(1,2) = -cofactor_4x4<MatrixType,2,1>(matrix);
00222     result.coeffRef(2,2) =  cofactor_4x4<MatrixType,2,2>(matrix);
00223     result.coeffRef(3,2) = -cofactor_4x4<MatrixType,2,3>(matrix);
00224     result.coeffRef(0,1) = -cofactor_4x4<MatrixType,1,0>(matrix);
00225     result.coeffRef(1,1) =  cofactor_4x4<MatrixType,1,1>(matrix);
00226     result.coeffRef(2,1) = -cofactor_4x4<MatrixType,1,2>(matrix);
00227     result.coeffRef(3,1) =  cofactor_4x4<MatrixType,1,3>(matrix);
00228     result.coeffRef(0,3) = -cofactor_4x4<MatrixType,3,0>(matrix);
00229     result.coeffRef(1,3) =  cofactor_4x4<MatrixType,3,1>(matrix);
00230     result.coeffRef(2,3) = -cofactor_4x4<MatrixType,3,2>(matrix);
00231     result.coeffRef(3,3) =  cofactor_4x4<MatrixType,3,3>(matrix);
00232     result /= (matrix.col(0).cwiseProduct(result.row(0).transpose())).sum();
00233   }
00234 };
00235 
00236 template<typename MatrixType, typename ResultType>
00237 struct compute_inverse<MatrixType, ResultType, 4>
00238  : compute_inverse_size4<Architecture::Target, typename MatrixType::Scalar,
00239                             MatrixType, ResultType>
00240 {
00241 };
00242 
00243 template<typename MatrixType, typename ResultType>
00244 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 4>
00245 {
00246   static inline void run(
00247     const MatrixType& matrix,
00248     const typename MatrixType::RealScalar& absDeterminantThreshold,
00249     ResultType& inverse,
00250     typename ResultType::Scalar& determinant,
00251     bool& invertible
00252   )
00253   {
00254     determinant = matrix.determinant();
00255     invertible = abs(determinant) > absDeterminantThreshold;
00256     if(invertible) compute_inverse<MatrixType, ResultType>::run(matrix, inverse);
00257   }
00258 };
00259 
00260 /*************************
00261 *** MatrixBase methods ***
00262 *************************/
00263 
00264 template<typename MatrixType>
00265 struct traits<inverse_impl<MatrixType> >
00266 {
00267   typedef typename MatrixType::PlainObject ReturnType;
00268 };
00269 
00270 template<typename MatrixType>
00271 struct inverse_impl : public ReturnByValue<inverse_impl<MatrixType> >
00272 {
00273   typedef typename MatrixType::Index Index;
00274   typedef typename internal::eval<MatrixType>::type MatrixTypeNested;
00275   typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
00276   MatrixTypeNested m_matrix;
00277 
00278   inverse_impl(const MatrixType& matrix)
00279     : m_matrix(matrix)
00280   {}
00281 
00282   inline Index rows() const { return m_matrix.rows(); }
00283   inline Index cols() const { return m_matrix.cols(); }
00284 
00285   template<typename Dest> inline void evalTo(Dest& dst) const
00286   {
00287     const int Size = EIGEN_PLAIN_ENUM_MIN(MatrixType::ColsAtCompileTime,Dest::ColsAtCompileTime);
00288     EIGEN_ONLY_USED_FOR_DEBUG(Size);
00289     eigen_assert(( (Size<=1) || (Size>4) || (extract_data(m_matrix)!=extract_data(dst)))
00290               && "Aliasing problem detected in inverse(), you need to do inverse().eval() here.");
00291 
00292     compute_inverse<MatrixTypeNestedCleaned, Dest>::run(m_matrix, dst);
00293   }
00294 };
00295 
00296 } // end namespace internal
00297 
00315 template<typename Derived>
00316 inline const internal::inverse_impl<Derived> MatrixBase<Derived>::inverse() const
00317 {
00318   EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsInteger,THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES)
00319   eigen_assert(rows() == cols());
00320   return internal::inverse_impl<Derived>(derived());
00321 }
00322 
00341 template<typename Derived>
00342 template<typename ResultType>
00343 inline void MatrixBase<Derived>::computeInverseAndDetWithCheck(
00344     ResultType& inverse,
00345     typename ResultType::Scalar& determinant,
00346     bool& invertible,
00347     const RealScalar& absDeterminantThreshold
00348   ) const
00349 {
00350   // i'd love to put some static assertions there, but SFINAE means that they have no effect...
00351   eigen_assert(rows() == cols());
00352   // for 2x2, it's worth giving a chance to avoid evaluating.
00353   // for larger sizes, evaluating has negligible cost and limits code size.
00354   typedef typename internal::conditional<
00355     RowsAtCompileTime == 2,
00356     typename internal::remove_all<typename internal::nested<Derived, 2>::type>::type,
00357     PlainObject
00358   >::type MatrixType;
00359   internal::compute_inverse_and_det_with_check<MatrixType, ResultType>::run
00360     (derived(), absDeterminantThreshold, inverse, determinant, invertible);
00361 }
00362 
00380 template<typename Derived>
00381 template<typename ResultType>
00382 inline void MatrixBase<Derived>::computeInverseWithCheck(
00383     ResultType& inverse,
00384     bool& invertible,
00385     const RealScalar& absDeterminantThreshold
00386   ) const
00387 {
00388   RealScalar determinant;
00389   // i'd love to put some static assertions there, but SFINAE means that they have no effect...
00390   eigen_assert(rows() == cols());
00391   computeInverseAndDetWithCheck(inverse,determinant,invertible,absDeterminantThreshold);
00392 }
00393 
00394 } // end namespace Eigen
00395 
00396 #endif // EIGEN_INVERSE_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:49