Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
InverseImpl.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_INVERSE_IMPL_H
12 #define EIGEN_INVERSE_IMPL_H
13 
14 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18 namespace internal {
19 
20 /**********************************
21 *** General case implementation ***
22 **********************************/
23 
24 template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
25 struct compute_inverse
26 {
27  EIGEN_DEVICE_FUNC
28  static inline void run(const MatrixType& matrix, ResultType& result)
29  {
30  result = matrix.partialPivLu().inverse();
31  }
32 };
33 
34 template<typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
35 struct compute_inverse_and_det_with_check { /* nothing! general case not supported. */ };
36 
37 /****************************
38 *** Size 1 implementation ***
39 ****************************/
40 
41 template<typename MatrixType, typename ResultType>
42 struct compute_inverse<MatrixType, ResultType, 1>
43 {
44  EIGEN_DEVICE_FUNC
45  static inline void run(const MatrixType& matrix, ResultType& result)
46  {
47  typedef typename MatrixType::Scalar Scalar;
48  internal::evaluator<MatrixType> matrixEval(matrix);
49  result.coeffRef(0,0) = Scalar(1) / matrixEval.coeff(0,0);
50  }
51 };
52 
53 template<typename MatrixType, typename ResultType>
54 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 1>
55 {
56  EIGEN_DEVICE_FUNC
57  static inline void run(
58  const MatrixType& matrix,
59  const typename MatrixType::RealScalar& absDeterminantThreshold,
60  ResultType& result,
61  typename ResultType::Scalar& determinant,
62  bool& invertible
63  )
64  {
65  using std::abs;
66  determinant = matrix.coeff(0,0);
67  invertible = abs(determinant) > absDeterminantThreshold;
68  if(invertible) result.coeffRef(0,0) = typename ResultType::Scalar(1) / determinant;
69  }
70 };
71 
72 /****************************
73 *** Size 2 implementation ***
74 ****************************/
75 
76 template<typename MatrixType, typename ResultType>
77 EIGEN_DEVICE_FUNC
78 inline void compute_inverse_size2_helper(
79  const MatrixType& matrix, const typename ResultType::Scalar& invdet,
80  ResultType& result)
81 {
82  typename ResultType::Scalar temp = matrix.coeff(0,0);
83  result.coeffRef(0,0) = matrix.coeff(1,1) * invdet;
84  result.coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
85  result.coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
86  result.coeffRef(1,1) = temp * invdet;
87 }
88 
89 template<typename MatrixType, typename ResultType>
90 struct compute_inverse<MatrixType, ResultType, 2>
91 {
92  EIGEN_DEVICE_FUNC
93  static inline void run(const MatrixType& matrix, ResultType& result)
94  {
95  typedef typename ResultType::Scalar Scalar;
96  const Scalar invdet = typename MatrixType::Scalar(1) / matrix.determinant();
97  compute_inverse_size2_helper(matrix, invdet, result);
98  }
99 };
100 
101 template<typename MatrixType, typename ResultType>
102 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 2>
103 {
104  EIGEN_DEVICE_FUNC
105  static inline void run(
106  const MatrixType& matrix,
107  const typename MatrixType::RealScalar& absDeterminantThreshold,
108  ResultType& inverse,
109  typename ResultType::Scalar& determinant,
110  bool& invertible
111  )
112  {
113  using std::abs;
114  typedef typename ResultType::Scalar Scalar;
115  determinant = matrix.determinant();
116  invertible = abs(determinant) > absDeterminantThreshold;
117  if(!invertible) return;
118  const Scalar invdet = Scalar(1) / determinant;
119  compute_inverse_size2_helper(matrix, invdet, inverse);
120  }
121 };
122 
123 /****************************
124 *** Size 3 implementation ***
125 ****************************/
126 
127 template<typename MatrixType, int i, int j>
128 EIGEN_DEVICE_FUNC
129 inline typename MatrixType::Scalar cofactor_3x3(const MatrixType& m)
130 {
131  enum {
132  i1 = (i+1) % 3,
133  i2 = (i+2) % 3,
134  j1 = (j+1) % 3,
135  j2 = (j+2) % 3
136  };
137  return m.coeff(i1, j1) * m.coeff(i2, j2)
138  - m.coeff(i1, j2) * m.coeff(i2, j1);
139 }
140 
141 template<typename MatrixType, typename ResultType>
142 EIGEN_DEVICE_FUNC
143 inline void compute_inverse_size3_helper(
144  const MatrixType& matrix,
145  const typename ResultType::Scalar& invdet,
146  const Matrix<typename ResultType::Scalar,3,1>& cofactors_col0,
147  ResultType& result)
148 {
149  // Compute cofactors in a way that avoids aliasing issues.
150  typedef typename ResultType::Scalar Scalar;
151  const Scalar c01 = cofactor_3x3<MatrixType,0,1>(matrix) * invdet;
152  const Scalar c11 = cofactor_3x3<MatrixType,1,1>(matrix) * invdet;
153  const Scalar c02 = cofactor_3x3<MatrixType,0,2>(matrix) * invdet;
154  result.coeffRef(1,2) = cofactor_3x3<MatrixType,2,1>(matrix) * invdet;
155  result.coeffRef(2,1) = cofactor_3x3<MatrixType,1,2>(matrix) * invdet;
156  result.coeffRef(2,2) = cofactor_3x3<MatrixType,2,2>(matrix) * invdet;
157  result.coeffRef(1,0) = c01;
158  result.coeffRef(1,1) = c11;
159  result.coeffRef(2,0) = c02;
160  result.row(0) = cofactors_col0 * invdet;
161 }
162 
163 template<typename MatrixType, typename ResultType>
164 struct compute_inverse<MatrixType, ResultType, 3>
165 {
166  EIGEN_DEVICE_FUNC
167  static inline void run(const MatrixType& matrix, ResultType& result)
168  {
169  typedef typename ResultType::Scalar Scalar;
170  Matrix<typename MatrixType::Scalar,3,1> cofactors_col0;
171  cofactors_col0.coeffRef(0) = cofactor_3x3<MatrixType,0,0>(matrix);
172  cofactors_col0.coeffRef(1) = cofactor_3x3<MatrixType,1,0>(matrix);
173  cofactors_col0.coeffRef(2) = cofactor_3x3<MatrixType,2,0>(matrix);
174  const Scalar det = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
175  const Scalar invdet = Scalar(1) / det;
176  compute_inverse_size3_helper(matrix, invdet, cofactors_col0, result);
177  }
178 };
179 
180 template<typename MatrixType, typename ResultType>
181 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 3>
182 {
183  EIGEN_DEVICE_FUNC
184  static inline void run(
185  const MatrixType& matrix,
186  const typename MatrixType::RealScalar& absDeterminantThreshold,
187  ResultType& inverse,
188  typename ResultType::Scalar& determinant,
189  bool& invertible
190  )
191  {
192  typedef typename ResultType::Scalar Scalar;
193  Matrix<Scalar,3,1> cofactors_col0;
194  cofactors_col0.coeffRef(0) = cofactor_3x3<MatrixType,0,0>(matrix);
195  cofactors_col0.coeffRef(1) = cofactor_3x3<MatrixType,1,0>(matrix);
196  cofactors_col0.coeffRef(2) = cofactor_3x3<MatrixType,2,0>(matrix);
197  determinant = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
198  invertible = Eigen::numext::abs(determinant) > absDeterminantThreshold;
199  if(!invertible) return;
200  const Scalar invdet = Scalar(1) / determinant;
201  compute_inverse_size3_helper(matrix, invdet, cofactors_col0, inverse);
202  }
203 };
204 
205 /****************************
206 *** Size 4 implementation ***
207 ****************************/
208 
209 template<typename Derived>
210 EIGEN_DEVICE_FUNC
211 inline const typename Derived::Scalar general_det3_helper
212 (const MatrixBase<Derived>& matrix, int i1, int i2, int i3, int j1, int j2, int j3)
213 {
214  return matrix.coeff(i1,j1)
215  * (matrix.coeff(i2,j2) * matrix.coeff(i3,j3) - matrix.coeff(i2,j3) * matrix.coeff(i3,j2));
216 }
217 
218 template<typename MatrixType, int i, int j>
219 EIGEN_DEVICE_FUNC
220 inline typename MatrixType::Scalar cofactor_4x4(const MatrixType& matrix)
221 {
222  enum {
223  i1 = (i+1) % 4,
224  i2 = (i+2) % 4,
225  i3 = (i+3) % 4,
226  j1 = (j+1) % 4,
227  j2 = (j+2) % 4,
228  j3 = (j+3) % 4
229  };
230  return general_det3_helper(matrix, i1, i2, i3, j1, j2, j3)
231  + general_det3_helper(matrix, i2, i3, i1, j1, j2, j3)
232  + general_det3_helper(matrix, i3, i1, i2, j1, j2, j3);
233 }
234 
235 template<int Arch, typename Scalar, typename MatrixType, typename ResultType>
236 struct compute_inverse_size4
237 {
238  EIGEN_DEVICE_FUNC
239  static void run(const MatrixType& matrix, ResultType& result)
240  {
241  result.coeffRef(0,0) = cofactor_4x4<MatrixType,0,0>(matrix);
242  result.coeffRef(1,0) = -cofactor_4x4<MatrixType,0,1>(matrix);
243  result.coeffRef(2,0) = cofactor_4x4<MatrixType,0,2>(matrix);
244  result.coeffRef(3,0) = -cofactor_4x4<MatrixType,0,3>(matrix);
245  result.coeffRef(0,2) = cofactor_4x4<MatrixType,2,0>(matrix);
246  result.coeffRef(1,2) = -cofactor_4x4<MatrixType,2,1>(matrix);
247  result.coeffRef(2,2) = cofactor_4x4<MatrixType,2,2>(matrix);
248  result.coeffRef(3,2) = -cofactor_4x4<MatrixType,2,3>(matrix);
249  result.coeffRef(0,1) = -cofactor_4x4<MatrixType,1,0>(matrix);
250  result.coeffRef(1,1) = cofactor_4x4<MatrixType,1,1>(matrix);
251  result.coeffRef(2,1) = -cofactor_4x4<MatrixType,1,2>(matrix);
252  result.coeffRef(3,1) = cofactor_4x4<MatrixType,1,3>(matrix);
253  result.coeffRef(0,3) = -cofactor_4x4<MatrixType,3,0>(matrix);
254  result.coeffRef(1,3) = cofactor_4x4<MatrixType,3,1>(matrix);
255  result.coeffRef(2,3) = -cofactor_4x4<MatrixType,3,2>(matrix);
256  result.coeffRef(3,3) = cofactor_4x4<MatrixType,3,3>(matrix);
257  result /= (matrix.col(0).cwiseProduct(result.row(0).transpose())).sum();
258  }
259 };
260 
261 template<typename MatrixType, typename ResultType>
262 struct compute_inverse<MatrixType, ResultType, 4>
263  : compute_inverse_size4<Architecture::Target, typename MatrixType::Scalar,
264  MatrixType, ResultType>
265 {
266 };
267 
268 template<typename MatrixType, typename ResultType>
269 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 4>
270 {
271  EIGEN_DEVICE_FUNC
272  static inline void run(
273  const MatrixType& matrix,
274  const typename MatrixType::RealScalar& absDeterminantThreshold,
275  ResultType& inverse,
276  typename ResultType::Scalar& determinant,
277  bool& invertible
278  )
279  {
280  using std::abs;
281  determinant = matrix.determinant();
282  invertible = abs(determinant) > absDeterminantThreshold;
283  if(invertible && extract_data(matrix) != extract_data(inverse)) {
284  compute_inverse<MatrixType, ResultType>::run(matrix, inverse);
285  }
286  else if(invertible) {
287  MatrixType matrix_t = matrix;
288  compute_inverse<MatrixType, ResultType>::run(matrix_t, inverse);
289  }
290  }
291 };
292 
293 /*************************
294 *** MatrixBase methods ***
295 *************************/
296 
297 } // end namespace internal
298 
299 namespace internal {
300 
301 // Specialization for "dense = dense_xpr.inverse()"
302 template<typename DstXprType, typename XprType>
303 struct Assignment<DstXprType, Inverse<XprType>, internal::assign_op<typename DstXprType::Scalar,typename XprType::Scalar>, Dense2Dense>
304 {
305  typedef Inverse<XprType> SrcXprType;
306  EIGEN_DEVICE_FUNC
307  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename XprType::Scalar> &)
308  {
309  Index dstRows = src.rows();
310  Index dstCols = src.cols();
311  if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
312  dst.resize(dstRows, dstCols);
313 
314  const int Size = plain_enum_min(XprType::ColsAtCompileTime, DstXprType::ColsAtCompileTime);
315  EIGEN_ONLY_USED_FOR_DEBUG(Size);
316  eigen_assert(( (Size<=1) || (Size>4) || (extract_data(src.nestedExpression())!=extract_data(dst)))
317  && "Aliasing problem detected in inverse(), you need to do inverse().eval() here.");
318 
319  typedef typename internal::nested_eval<XprType,XprType::ColsAtCompileTime>::type ActualXprType;
320  typedef internal::remove_all_t<ActualXprType> ActualXprTypeCleanded;
321 
322  ActualXprType actual_xpr(src.nestedExpression());
323 
324  compute_inverse<ActualXprTypeCleanded, DstXprType>::run(actual_xpr, dst);
325  }
326 };
327 
328 
329 } // end namespace internal
330 
348 template<typename Derived>
349 EIGEN_DEVICE_FUNC
351 {
352  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsInteger,THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES)
353  eigen_assert(rows() == cols());
354  return Inverse<Derived>(derived());
355 }
356 
377 template<typename Derived>
378 template<typename ResultType>
380  ResultType& inverse,
381  typename ResultType::Scalar& determinant,
382  bool& invertible,
383  const RealScalar& absDeterminantThreshold
384  ) const
385 {
386  // i'd love to put some static assertions there, but SFINAE means that they have no effect...
387  eigen_assert(rows() == cols());
388  // for 2x2, it's worth giving a chance to avoid evaluating.
389  // for larger sizes, evaluating has negligible cost and limits code size.
390  typedef std::conditional_t<
391  RowsAtCompileTime == 2,
392  internal::remove_all_t<typename internal::nested_eval<Derived, 2>::type>,
394  > MatrixType;
395  internal::compute_inverse_and_det_with_check<MatrixType, ResultType>::run
396  (derived(), absDeterminantThreshold, inverse, determinant, invertible);
397 }
398 
418 template<typename Derived>
419 template<typename ResultType>
421  ResultType& inverse,
422  bool& invertible,
423  const RealScalar& absDeterminantThreshold
424  ) const
425 {
426  Scalar determinant;
427  // i'd love to put some static assertions there, but SFINAE means that they have no effect...
428  eigen_assert(rows() == cols());
429  computeInverseAndDetWithCheck(inverse,determinant,invertible,absDeterminantThreshold);
430 }
431 
432 } // end namespace Eigen
433 
434 #endif // EIGEN_INVERSE_IMPL_H
std::conditional_t< internal::is_same< typename internal::traits< Derived >::XprKind, MatrixXpr >::value, PlainMatrix, PlainArray > PlainObject
The plain matrix or array type corresponding to this expression.
Definition: DenseBase.h:204
internal::traits< Derived >::Scalar Scalar
Definition: DenseBase.h:61
Expression of the inverse of another expression.
Definition: Inverse.h:46
void computeInverseWithCheck(ResultType &inverse, bool &invertible, const RealScalar &absDeterminantThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: InverseImpl.h:420
const Inverse< Derived > inverse() const
Definition: InverseImpl.h:350
void computeInverseAndDetWithCheck(ResultType &inverse, typename ResultType::Scalar &determinant, bool &invertible, const RealScalar &absDeterminantThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: InverseImpl.h:379
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_inverse_op< typename Derived::Scalar >, const Derived > inverse(const Eigen::ArrayBase< Derived > &x)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231