Michael Ernst Peter / Eigen

Eigen libary for mbed

src/LU/Inverse.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 2008-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
ykuroda 0:13a5d365ba16 5 //
ykuroda 0:13a5d365ba16 6 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 7 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 9
ykuroda 0:13a5d365ba16 10 #ifndef EIGEN_INVERSE_H
ykuroda 0:13a5d365ba16 11 #define EIGEN_INVERSE_H
ykuroda 0:13a5d365ba16 12
ykuroda 0:13a5d365ba16 13 namespace Eigen {
ykuroda 0:13a5d365ba16 14
ykuroda 0:13a5d365ba16 15 namespace internal {
ykuroda 0:13a5d365ba16 16
ykuroda 0:13a5d365ba16 17 /**********************************
ykuroda 0:13a5d365ba16 18 *** General case implementation ***
ykuroda 0:13a5d365ba16 19 **********************************/
ykuroda 0:13a5d365ba16 20
ykuroda 0:13a5d365ba16 21 template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
ykuroda 0:13a5d365ba16 22 struct compute_inverse
ykuroda 0:13a5d365ba16 23 {
ykuroda 0:13a5d365ba16 24 static inline void run(const MatrixType& matrix, ResultType& result)
ykuroda 0:13a5d365ba16 25 {
ykuroda 0:13a5d365ba16 26 result = matrix.partialPivLu().inverse();
ykuroda 0:13a5d365ba16 27 }
ykuroda 0:13a5d365ba16 28 };
ykuroda 0:13a5d365ba16 29
ykuroda 0:13a5d365ba16 30 template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
ykuroda 0:13a5d365ba16 31 struct compute_inverse_and_det_with_check { /* nothing! general case not supported. */ };
ykuroda 0:13a5d365ba16 32
ykuroda 0:13a5d365ba16 33 /****************************
ykuroda 0:13a5d365ba16 34 *** Size 1 implementation ***
ykuroda 0:13a5d365ba16 35 ****************************/
ykuroda 0:13a5d365ba16 36
ykuroda 0:13a5d365ba16 37 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 38 struct compute_inverse<MatrixType, ResultType, 1>
ykuroda 0:13a5d365ba16 39 {
ykuroda 0:13a5d365ba16 40 static inline void run(const MatrixType& matrix, ResultType& result)
ykuroda 0:13a5d365ba16 41 {
ykuroda 0:13a5d365ba16 42 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 43 result.coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
ykuroda 0:13a5d365ba16 44 }
ykuroda 0:13a5d365ba16 45 };
ykuroda 0:13a5d365ba16 46
ykuroda 0:13a5d365ba16 47 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 48 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 1>
ykuroda 0:13a5d365ba16 49 {
ykuroda 0:13a5d365ba16 50 static inline void run(
ykuroda 0:13a5d365ba16 51 const MatrixType& matrix,
ykuroda 0:13a5d365ba16 52 const typename MatrixType::RealScalar& absDeterminantThreshold,
ykuroda 0:13a5d365ba16 53 ResultType& result,
ykuroda 0:13a5d365ba16 54 typename ResultType::Scalar& determinant,
ykuroda 0:13a5d365ba16 55 bool& invertible
ykuroda 0:13a5d365ba16 56 )
ykuroda 0:13a5d365ba16 57 {
ykuroda 0:13a5d365ba16 58 using std::abs;
ykuroda 0:13a5d365ba16 59 determinant = matrix.coeff(0,0);
ykuroda 0:13a5d365ba16 60 invertible = abs(determinant) > absDeterminantThreshold;
ykuroda 0:13a5d365ba16 61 if(invertible) result.coeffRef(0,0) = typename ResultType::Scalar(1) / determinant;
ykuroda 0:13a5d365ba16 62 }
ykuroda 0:13a5d365ba16 63 };
ykuroda 0:13a5d365ba16 64
ykuroda 0:13a5d365ba16 65 /****************************
ykuroda 0:13a5d365ba16 66 *** Size 2 implementation ***
ykuroda 0:13a5d365ba16 67 ****************************/
ykuroda 0:13a5d365ba16 68
ykuroda 0:13a5d365ba16 69 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 70 inline void compute_inverse_size2_helper(
ykuroda 0:13a5d365ba16 71 const MatrixType& matrix, const typename ResultType::Scalar& invdet,
ykuroda 0:13a5d365ba16 72 ResultType& result)
ykuroda 0:13a5d365ba16 73 {
ykuroda 0:13a5d365ba16 74 result.coeffRef(0,0) = matrix.coeff(1,1) * invdet;
ykuroda 0:13a5d365ba16 75 result.coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
ykuroda 0:13a5d365ba16 76 result.coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
ykuroda 0:13a5d365ba16 77 result.coeffRef(1,1) = matrix.coeff(0,0) * invdet;
ykuroda 0:13a5d365ba16 78 }
ykuroda 0:13a5d365ba16 79
ykuroda 0:13a5d365ba16 80 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 81 struct compute_inverse<MatrixType, ResultType, 2>
ykuroda 0:13a5d365ba16 82 {
ykuroda 0:13a5d365ba16 83 static inline void run(const MatrixType& matrix, ResultType& result)
ykuroda 0:13a5d365ba16 84 {
ykuroda 0:13a5d365ba16 85 typedef typename ResultType::Scalar Scalar;
ykuroda 0:13a5d365ba16 86 const Scalar invdet = typename MatrixType::Scalar(1) / matrix.determinant();
ykuroda 0:13a5d365ba16 87 compute_inverse_size2_helper(matrix, invdet, result);
ykuroda 0:13a5d365ba16 88 }
ykuroda 0:13a5d365ba16 89 };
ykuroda 0:13a5d365ba16 90
ykuroda 0:13a5d365ba16 91 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 92 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 2>
ykuroda 0:13a5d365ba16 93 {
ykuroda 0:13a5d365ba16 94 static inline void run(
ykuroda 0:13a5d365ba16 95 const MatrixType& matrix,
ykuroda 0:13a5d365ba16 96 const typename MatrixType::RealScalar& absDeterminantThreshold,
ykuroda 0:13a5d365ba16 97 ResultType& inverse,
ykuroda 0:13a5d365ba16 98 typename ResultType::Scalar& determinant,
ykuroda 0:13a5d365ba16 99 bool& invertible
ykuroda 0:13a5d365ba16 100 )
ykuroda 0:13a5d365ba16 101 {
ykuroda 0:13a5d365ba16 102 using std::abs;
ykuroda 0:13a5d365ba16 103 typedef typename ResultType::Scalar Scalar;
ykuroda 0:13a5d365ba16 104 determinant = matrix.determinant();
ykuroda 0:13a5d365ba16 105 invertible = abs(determinant) > absDeterminantThreshold;
ykuroda 0:13a5d365ba16 106 if(!invertible) return;
ykuroda 0:13a5d365ba16 107 const Scalar invdet = Scalar(1) / determinant;
ykuroda 0:13a5d365ba16 108 compute_inverse_size2_helper(matrix, invdet, inverse);
ykuroda 0:13a5d365ba16 109 }
ykuroda 0:13a5d365ba16 110 };
ykuroda 0:13a5d365ba16 111
ykuroda 0:13a5d365ba16 112 /****************************
ykuroda 0:13a5d365ba16 113 *** Size 3 implementation ***
ykuroda 0:13a5d365ba16 114 ****************************/
ykuroda 0:13a5d365ba16 115
ykuroda 0:13a5d365ba16 116 template<typename MatrixType, int i, int j>
ykuroda 0:13a5d365ba16 117 inline typename MatrixType::Scalar cofactor_3x3(const MatrixType& m)
ykuroda 0:13a5d365ba16 118 {
ykuroda 0:13a5d365ba16 119 enum {
ykuroda 0:13a5d365ba16 120 i1 = (i+1) % 3,
ykuroda 0:13a5d365ba16 121 i2 = (i+2) % 3,
ykuroda 0:13a5d365ba16 122 j1 = (j+1) % 3,
ykuroda 0:13a5d365ba16 123 j2 = (j+2) % 3
ykuroda 0:13a5d365ba16 124 };
ykuroda 0:13a5d365ba16 125 return m.coeff(i1, j1) * m.coeff(i2, j2)
ykuroda 0:13a5d365ba16 126 - m.coeff(i1, j2) * m.coeff(i2, j1);
ykuroda 0:13a5d365ba16 127 }
ykuroda 0:13a5d365ba16 128
ykuroda 0:13a5d365ba16 129 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 130 inline void compute_inverse_size3_helper(
ykuroda 0:13a5d365ba16 131 const MatrixType& matrix,
ykuroda 0:13a5d365ba16 132 const typename ResultType::Scalar& invdet,
ykuroda 0:13a5d365ba16 133 const Matrix<typename ResultType::Scalar,3,1>& cofactors_col0,
ykuroda 0:13a5d365ba16 134 ResultType& result)
ykuroda 0:13a5d365ba16 135 {
ykuroda 0:13a5d365ba16 136 result.row(0) = cofactors_col0 * invdet;
ykuroda 0:13a5d365ba16 137 result.coeffRef(1,0) = cofactor_3x3<MatrixType,0,1>(matrix) * invdet;
ykuroda 0:13a5d365ba16 138 result.coeffRef(1,1) = cofactor_3x3<MatrixType,1,1>(matrix) * invdet;
ykuroda 0:13a5d365ba16 139 result.coeffRef(1,2) = cofactor_3x3<MatrixType,2,1>(matrix) * invdet;
ykuroda 0:13a5d365ba16 140 result.coeffRef(2,0) = cofactor_3x3<MatrixType,0,2>(matrix) * invdet;
ykuroda 0:13a5d365ba16 141 result.coeffRef(2,1) = cofactor_3x3<MatrixType,1,2>(matrix) * invdet;
ykuroda 0:13a5d365ba16 142 result.coeffRef(2,2) = cofactor_3x3<MatrixType,2,2>(matrix) * invdet;
ykuroda 0:13a5d365ba16 143 }
ykuroda 0:13a5d365ba16 144
ykuroda 0:13a5d365ba16 145 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 146 struct compute_inverse<MatrixType, ResultType, 3>
ykuroda 0:13a5d365ba16 147 {
ykuroda 0:13a5d365ba16 148 static inline void run(const MatrixType& matrix, ResultType& result)
ykuroda 0:13a5d365ba16 149 {
ykuroda 0:13a5d365ba16 150 typedef typename ResultType::Scalar Scalar;
ykuroda 0:13a5d365ba16 151 Matrix<typename MatrixType::Scalar,3,1> cofactors_col0;
ykuroda 0:13a5d365ba16 152 cofactors_col0.coeffRef(0) = cofactor_3x3<MatrixType,0,0>(matrix);
ykuroda 0:13a5d365ba16 153 cofactors_col0.coeffRef(1) = cofactor_3x3<MatrixType,1,0>(matrix);
ykuroda 0:13a5d365ba16 154 cofactors_col0.coeffRef(2) = cofactor_3x3<MatrixType,2,0>(matrix);
ykuroda 0:13a5d365ba16 155 const Scalar det = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
ykuroda 0:13a5d365ba16 156 const Scalar invdet = Scalar(1) / det;
ykuroda 0:13a5d365ba16 157 compute_inverse_size3_helper(matrix, invdet, cofactors_col0, result);
ykuroda 0:13a5d365ba16 158 }
ykuroda 0:13a5d365ba16 159 };
ykuroda 0:13a5d365ba16 160
ykuroda 0:13a5d365ba16 161 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 162 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 3>
ykuroda 0:13a5d365ba16 163 {
ykuroda 0:13a5d365ba16 164 static inline void run(
ykuroda 0:13a5d365ba16 165 const MatrixType& matrix,
ykuroda 0:13a5d365ba16 166 const typename MatrixType::RealScalar& absDeterminantThreshold,
ykuroda 0:13a5d365ba16 167 ResultType& inverse,
ykuroda 0:13a5d365ba16 168 typename ResultType::Scalar& determinant,
ykuroda 0:13a5d365ba16 169 bool& invertible
ykuroda 0:13a5d365ba16 170 )
ykuroda 0:13a5d365ba16 171 {
ykuroda 0:13a5d365ba16 172 using std::abs;
ykuroda 0:13a5d365ba16 173 typedef typename ResultType::Scalar Scalar;
ykuroda 0:13a5d365ba16 174 Matrix<Scalar,3,1> cofactors_col0;
ykuroda 0:13a5d365ba16 175 cofactors_col0.coeffRef(0) = cofactor_3x3<MatrixType,0,0>(matrix);
ykuroda 0:13a5d365ba16 176 cofactors_col0.coeffRef(1) = cofactor_3x3<MatrixType,1,0>(matrix);
ykuroda 0:13a5d365ba16 177 cofactors_col0.coeffRef(2) = cofactor_3x3<MatrixType,2,0>(matrix);
ykuroda 0:13a5d365ba16 178 determinant = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
ykuroda 0:13a5d365ba16 179 invertible = abs(determinant) > absDeterminantThreshold;
ykuroda 0:13a5d365ba16 180 if(!invertible) return;
ykuroda 0:13a5d365ba16 181 const Scalar invdet = Scalar(1) / determinant;
ykuroda 0:13a5d365ba16 182 compute_inverse_size3_helper(matrix, invdet, cofactors_col0, inverse);
ykuroda 0:13a5d365ba16 183 }
ykuroda 0:13a5d365ba16 184 };
ykuroda 0:13a5d365ba16 185
ykuroda 0:13a5d365ba16 186 /****************************
ykuroda 0:13a5d365ba16 187 *** Size 4 implementation ***
ykuroda 0:13a5d365ba16 188 ****************************/
ykuroda 0:13a5d365ba16 189
ykuroda 0:13a5d365ba16 190 template<typename Derived>
ykuroda 0:13a5d365ba16 191 inline const typename Derived::Scalar general_det3_helper
ykuroda 0:13a5d365ba16 192 (const MatrixBase<Derived>& matrix, int i1, int i2, int i3, int j1, int j2, int j3)
ykuroda 0:13a5d365ba16 193 {
ykuroda 0:13a5d365ba16 194 return matrix.coeff(i1,j1)
ykuroda 0:13a5d365ba16 195 * (matrix.coeff(i2,j2) * matrix.coeff(i3,j3) - matrix.coeff(i2,j3) * matrix.coeff(i3,j2));
ykuroda 0:13a5d365ba16 196 }
ykuroda 0:13a5d365ba16 197
ykuroda 0:13a5d365ba16 198 template<typename MatrixType, int i, int j>
ykuroda 0:13a5d365ba16 199 inline typename MatrixType::Scalar cofactor_4x4(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 200 {
ykuroda 0:13a5d365ba16 201 enum {
ykuroda 0:13a5d365ba16 202 i1 = (i+1) % 4,
ykuroda 0:13a5d365ba16 203 i2 = (i+2) % 4,
ykuroda 0:13a5d365ba16 204 i3 = (i+3) % 4,
ykuroda 0:13a5d365ba16 205 j1 = (j+1) % 4,
ykuroda 0:13a5d365ba16 206 j2 = (j+2) % 4,
ykuroda 0:13a5d365ba16 207 j3 = (j+3) % 4
ykuroda 0:13a5d365ba16 208 };
ykuroda 0:13a5d365ba16 209 return general_det3_helper(matrix, i1, i2, i3, j1, j2, j3)
ykuroda 0:13a5d365ba16 210 + general_det3_helper(matrix, i2, i3, i1, j1, j2, j3)
ykuroda 0:13a5d365ba16 211 + general_det3_helper(matrix, i3, i1, i2, j1, j2, j3);
ykuroda 0:13a5d365ba16 212 }
ykuroda 0:13a5d365ba16 213
ykuroda 0:13a5d365ba16 214 template<int Arch, typename Scalar, typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 215 struct compute_inverse_size4
ykuroda 0:13a5d365ba16 216 {
ykuroda 0:13a5d365ba16 217 static void run(const MatrixType& matrix, ResultType& result)
ykuroda 0:13a5d365ba16 218 {
ykuroda 0:13a5d365ba16 219 result.coeffRef(0,0) = cofactor_4x4<MatrixType,0,0>(matrix);
ykuroda 0:13a5d365ba16 220 result.coeffRef(1,0) = -cofactor_4x4<MatrixType,0,1>(matrix);
ykuroda 0:13a5d365ba16 221 result.coeffRef(2,0) = cofactor_4x4<MatrixType,0,2>(matrix);
ykuroda 0:13a5d365ba16 222 result.coeffRef(3,0) = -cofactor_4x4<MatrixType,0,3>(matrix);
ykuroda 0:13a5d365ba16 223 result.coeffRef(0,2) = cofactor_4x4<MatrixType,2,0>(matrix);
ykuroda 0:13a5d365ba16 224 result.coeffRef(1,2) = -cofactor_4x4<MatrixType,2,1>(matrix);
ykuroda 0:13a5d365ba16 225 result.coeffRef(2,2) = cofactor_4x4<MatrixType,2,2>(matrix);
ykuroda 0:13a5d365ba16 226 result.coeffRef(3,2) = -cofactor_4x4<MatrixType,2,3>(matrix);
ykuroda 0:13a5d365ba16 227 result.coeffRef(0,1) = -cofactor_4x4<MatrixType,1,0>(matrix);
ykuroda 0:13a5d365ba16 228 result.coeffRef(1,1) = cofactor_4x4<MatrixType,1,1>(matrix);
ykuroda 0:13a5d365ba16 229 result.coeffRef(2,1) = -cofactor_4x4<MatrixType,1,2>(matrix);
ykuroda 0:13a5d365ba16 230 result.coeffRef(3,1) = cofactor_4x4<MatrixType,1,3>(matrix);
ykuroda 0:13a5d365ba16 231 result.coeffRef(0,3) = -cofactor_4x4<MatrixType,3,0>(matrix);
ykuroda 0:13a5d365ba16 232 result.coeffRef(1,3) = cofactor_4x4<MatrixType,3,1>(matrix);
ykuroda 0:13a5d365ba16 233 result.coeffRef(2,3) = -cofactor_4x4<MatrixType,3,2>(matrix);
ykuroda 0:13a5d365ba16 234 result.coeffRef(3,3) = cofactor_4x4<MatrixType,3,3>(matrix);
ykuroda 0:13a5d365ba16 235 result /= (matrix.col(0).cwiseProduct(result.row(0).transpose())).sum();
ykuroda 0:13a5d365ba16 236 }
ykuroda 0:13a5d365ba16 237 };
ykuroda 0:13a5d365ba16 238
ykuroda 0:13a5d365ba16 239 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 240 struct compute_inverse<MatrixType, ResultType, 4>
ykuroda 0:13a5d365ba16 241 : compute_inverse_size4<Architecture::Target, typename MatrixType::Scalar,
ykuroda 0:13a5d365ba16 242 MatrixType, ResultType>
ykuroda 0:13a5d365ba16 243 {
ykuroda 0:13a5d365ba16 244 };
ykuroda 0:13a5d365ba16 245
ykuroda 0:13a5d365ba16 246 template<typename MatrixType, typename ResultType>
ykuroda 0:13a5d365ba16 247 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 4>
ykuroda 0:13a5d365ba16 248 {
ykuroda 0:13a5d365ba16 249 static inline void run(
ykuroda 0:13a5d365ba16 250 const MatrixType& matrix,
ykuroda 0:13a5d365ba16 251 const typename MatrixType::RealScalar& absDeterminantThreshold,
ykuroda 0:13a5d365ba16 252 ResultType& inverse,
ykuroda 0:13a5d365ba16 253 typename ResultType::Scalar& determinant,
ykuroda 0:13a5d365ba16 254 bool& invertible
ykuroda 0:13a5d365ba16 255 )
ykuroda 0:13a5d365ba16 256 {
ykuroda 0:13a5d365ba16 257 using std::abs;
ykuroda 0:13a5d365ba16 258 determinant = matrix.determinant();
ykuroda 0:13a5d365ba16 259 invertible = abs(determinant) > absDeterminantThreshold;
ykuroda 0:13a5d365ba16 260 if(invertible) compute_inverse<MatrixType, ResultType>::run(matrix, inverse);
ykuroda 0:13a5d365ba16 261 }
ykuroda 0:13a5d365ba16 262 };
ykuroda 0:13a5d365ba16 263
ykuroda 0:13a5d365ba16 264 /*************************
ykuroda 0:13a5d365ba16 265 *** MatrixBase methods ***
ykuroda 0:13a5d365ba16 266 *************************/
ykuroda 0:13a5d365ba16 267
ykuroda 0:13a5d365ba16 268 template<typename MatrixType>
ykuroda 0:13a5d365ba16 269 struct traits<inverse_impl<MatrixType> >
ykuroda 0:13a5d365ba16 270 {
ykuroda 0:13a5d365ba16 271 typedef typename MatrixType::PlainObject ReturnType;
ykuroda 0:13a5d365ba16 272 };
ykuroda 0:13a5d365ba16 273
ykuroda 0:13a5d365ba16 274 template<typename MatrixType>
ykuroda 0:13a5d365ba16 275 struct inverse_impl : public ReturnByValue<inverse_impl<MatrixType> >
ykuroda 0:13a5d365ba16 276 {
ykuroda 0:13a5d365ba16 277 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 278 typedef typename internal::eval<MatrixType>::type MatrixTypeNested;
ykuroda 0:13a5d365ba16 279 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
ykuroda 0:13a5d365ba16 280 MatrixTypeNested m_matrix;
ykuroda 0:13a5d365ba16 281
ykuroda 0:13a5d365ba16 282 inverse_impl(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 283 : m_matrix(matrix)
ykuroda 0:13a5d365ba16 284 {}
ykuroda 0:13a5d365ba16 285
ykuroda 0:13a5d365ba16 286 inline Index rows() const { return m_matrix.rows(); }
ykuroda 0:13a5d365ba16 287 inline Index cols() const { return m_matrix.cols(); }
ykuroda 0:13a5d365ba16 288
ykuroda 0:13a5d365ba16 289 template<typename Dest> inline void evalTo(Dest& dst) const
ykuroda 0:13a5d365ba16 290 {
ykuroda 0:13a5d365ba16 291 const int Size = EIGEN_PLAIN_ENUM_MIN(MatrixType::ColsAtCompileTime,Dest::ColsAtCompileTime);
ykuroda 0:13a5d365ba16 292 EIGEN_ONLY_USED_FOR_DEBUG(Size);
ykuroda 0:13a5d365ba16 293 eigen_assert(( (Size<=1) || (Size>4) || (extract_data(m_matrix)!=extract_data(dst)))
ykuroda 0:13a5d365ba16 294 && "Aliasing problem detected in inverse(), you need to do inverse().eval() here.");
ykuroda 0:13a5d365ba16 295
ykuroda 0:13a5d365ba16 296 compute_inverse<MatrixTypeNestedCleaned, Dest>::run(m_matrix, dst);
ykuroda 0:13a5d365ba16 297 }
ykuroda 0:13a5d365ba16 298 };
ykuroda 0:13a5d365ba16 299
ykuroda 0:13a5d365ba16 300 } // end namespace internal
ykuroda 0:13a5d365ba16 301
ykuroda 0:13a5d365ba16 302 /** \lu_module
ykuroda 0:13a5d365ba16 303 *
ykuroda 0:13a5d365ba16 304 * \returns the matrix inverse of this matrix.
ykuroda 0:13a5d365ba16 305 *
ykuroda 0:13a5d365ba16 306 * For small fixed sizes up to 4x4, this method uses cofactors.
ykuroda 0:13a5d365ba16 307 * In the general case, this method uses class PartialPivLU.
ykuroda 0:13a5d365ba16 308 *
ykuroda 0:13a5d365ba16 309 * \note This matrix must be invertible, otherwise the result is undefined. If you need an
ykuroda 0:13a5d365ba16 310 * invertibility check, do the following:
ykuroda 0:13a5d365ba16 311 * \li for fixed sizes up to 4x4, use computeInverseAndDetWithCheck().
ykuroda 0:13a5d365ba16 312 * \li for the general case, use class FullPivLU.
ykuroda 0:13a5d365ba16 313 *
ykuroda 0:13a5d365ba16 314 * Example: \include MatrixBase_inverse.cpp
ykuroda 0:13a5d365ba16 315 * Output: \verbinclude MatrixBase_inverse.out
ykuroda 0:13a5d365ba16 316 *
ykuroda 0:13a5d365ba16 317 * \sa computeInverseAndDetWithCheck()
ykuroda 0:13a5d365ba16 318 */
ykuroda 0:13a5d365ba16 319 template<typename Derived>
ykuroda 0:13a5d365ba16 320 inline const internal::inverse_impl<Derived> MatrixBase<Derived>::inverse() const
ykuroda 0:13a5d365ba16 321 {
ykuroda 0:13a5d365ba16 322 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsInteger,THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES)
ykuroda 0:13a5d365ba16 323 eigen_assert(rows() == cols());
ykuroda 0:13a5d365ba16 324 return internal::inverse_impl<Derived>(derived());
ykuroda 0:13a5d365ba16 325 }
ykuroda 0:13a5d365ba16 326
ykuroda 0:13a5d365ba16 327 /** \lu_module
ykuroda 0:13a5d365ba16 328 *
ykuroda 0:13a5d365ba16 329 * Computation of matrix inverse and determinant, with invertibility check.
ykuroda 0:13a5d365ba16 330 *
ykuroda 0:13a5d365ba16 331 * This is only for fixed-size square matrices of size up to 4x4.
ykuroda 0:13a5d365ba16 332 *
ykuroda 0:13a5d365ba16 333 * \param inverse Reference to the matrix in which to store the inverse.
ykuroda 0:13a5d365ba16 334 * \param determinant Reference to the variable in which to store the determinant.
ykuroda 0:13a5d365ba16 335 * \param invertible Reference to the bool variable in which to store whether the matrix is invertible.
ykuroda 0:13a5d365ba16 336 * \param absDeterminantThreshold Optional parameter controlling the invertibility check.
ykuroda 0:13a5d365ba16 337 * The matrix will be declared invertible if the absolute value of its
ykuroda 0:13a5d365ba16 338 * determinant is greater than this threshold.
ykuroda 0:13a5d365ba16 339 *
ykuroda 0:13a5d365ba16 340 * Example: \include MatrixBase_computeInverseAndDetWithCheck.cpp
ykuroda 0:13a5d365ba16 341 * Output: \verbinclude MatrixBase_computeInverseAndDetWithCheck.out
ykuroda 0:13a5d365ba16 342 *
ykuroda 0:13a5d365ba16 343 * \sa inverse(), computeInverseWithCheck()
ykuroda 0:13a5d365ba16 344 */
ykuroda 0:13a5d365ba16 345 template<typename Derived>
ykuroda 0:13a5d365ba16 346 template<typename ResultType>
ykuroda 0:13a5d365ba16 347 inline void MatrixBase<Derived>::computeInverseAndDetWithCheck(
ykuroda 0:13a5d365ba16 348 ResultType& inverse,
ykuroda 0:13a5d365ba16 349 typename ResultType::Scalar& determinant,
ykuroda 0:13a5d365ba16 350 bool& invertible,
ykuroda 0:13a5d365ba16 351 const RealScalar& absDeterminantThreshold
ykuroda 0:13a5d365ba16 352 ) const
ykuroda 0:13a5d365ba16 353 {
ykuroda 0:13a5d365ba16 354 // i'd love to put some static assertions there, but SFINAE means that they have no effect...
ykuroda 0:13a5d365ba16 355 eigen_assert(rows() == cols());
ykuroda 0:13a5d365ba16 356 // for 2x2, it's worth giving a chance to avoid evaluating.
ykuroda 0:13a5d365ba16 357 // for larger sizes, evaluating has negligible cost and limits code size.
ykuroda 0:13a5d365ba16 358 typedef typename internal::conditional<
ykuroda 0:13a5d365ba16 359 RowsAtCompileTime == 2,
ykuroda 0:13a5d365ba16 360 typename internal::remove_all<typename internal::nested<Derived, 2>::type>::type,
ykuroda 0:13a5d365ba16 361 PlainObject
ykuroda 0:13a5d365ba16 362 >::type MatrixType;
ykuroda 0:13a5d365ba16 363 internal::compute_inverse_and_det_with_check<MatrixType, ResultType>::run
ykuroda 0:13a5d365ba16 364 (derived(), absDeterminantThreshold, inverse, determinant, invertible);
ykuroda 0:13a5d365ba16 365 }
ykuroda 0:13a5d365ba16 366
ykuroda 0:13a5d365ba16 367 /** \lu_module
ykuroda 0:13a5d365ba16 368 *
ykuroda 0:13a5d365ba16 369 * Computation of matrix inverse, with invertibility check.
ykuroda 0:13a5d365ba16 370 *
ykuroda 0:13a5d365ba16 371 * This is only for fixed-size square matrices of size up to 4x4.
ykuroda 0:13a5d365ba16 372 *
ykuroda 0:13a5d365ba16 373 * \param inverse Reference to the matrix in which to store the inverse.
ykuroda 0:13a5d365ba16 374 * \param invertible Reference to the bool variable in which to store whether the matrix is invertible.
ykuroda 0:13a5d365ba16 375 * \param absDeterminantThreshold Optional parameter controlling the invertibility check.
ykuroda 0:13a5d365ba16 376 * The matrix will be declared invertible if the absolute value of its
ykuroda 0:13a5d365ba16 377 * determinant is greater than this threshold.
ykuroda 0:13a5d365ba16 378 *
ykuroda 0:13a5d365ba16 379 * Example: \include MatrixBase_computeInverseWithCheck.cpp
ykuroda 0:13a5d365ba16 380 * Output: \verbinclude MatrixBase_computeInverseWithCheck.out
ykuroda 0:13a5d365ba16 381 *
ykuroda 0:13a5d365ba16 382 * \sa inverse(), computeInverseAndDetWithCheck()
ykuroda 0:13a5d365ba16 383 */
ykuroda 0:13a5d365ba16 384 template<typename Derived>
ykuroda 0:13a5d365ba16 385 template<typename ResultType>
ykuroda 0:13a5d365ba16 386 inline void MatrixBase<Derived>::computeInverseWithCheck(
ykuroda 0:13a5d365ba16 387 ResultType& inverse,
ykuroda 0:13a5d365ba16 388 bool& invertible,
ykuroda 0:13a5d365ba16 389 const RealScalar& absDeterminantThreshold
ykuroda 0:13a5d365ba16 390 ) const
ykuroda 0:13a5d365ba16 391 {
ykuroda 0:13a5d365ba16 392 RealScalar determinant;
ykuroda 0:13a5d365ba16 393 // i'd love to put some static assertions there, but SFINAE means that they have no effect...
ykuroda 0:13a5d365ba16 394 eigen_assert(rows() == cols());
ykuroda 0:13a5d365ba16 395 computeInverseAndDetWithCheck(inverse,determinant,invertible,absDeterminantThreshold);
ykuroda 0:13a5d365ba16 396 }
ykuroda 0:13a5d365ba16 397
ykuroda 0:13a5d365ba16 398 } // end namespace Eigen
ykuroda 0:13a5d365ba16 399
ykuroda 0:13a5d365ba16 400 #endif // EIGEN_INVERSE_H