Eigen  3.4.90 (git rev 5a9f66fb35d03a4da9ef8976e67a61b30aa16dcf)
 
Loading...
Searching...
No Matches
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 <[email protected]>
5// Copyright (C) 2014 Gael Guennebaud <[email protected]>
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// IWYU pragma: private
15#include "./InternalHeaderCheck.h"
16
17namespace Eigen {
18
19namespace internal {
20
21/**********************************
22*** General case implementation ***
23**********************************/
24
25template <typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
26struct compute_inverse {
27 EIGEN_DEVICE_FUNC static inline void run(const MatrixType& matrix, ResultType& result) {
28 result = matrix.partialPivLu().inverse();
29 }
30};
31
32template <typename MatrixType, typename ResultType, int Size = MatrixType::RowsAtCompileTime>
33struct compute_inverse_and_det_with_check { /* nothing! general case not supported. */
34};
35
36/****************************
37*** Size 1 implementation ***
38****************************/
39
40template <typename MatrixType, typename ResultType>
41struct compute_inverse<MatrixType, ResultType, 1> {
42 EIGEN_DEVICE_FUNC static inline void run(const MatrixType& matrix, ResultType& result) {
43 typedef typename MatrixType::Scalar Scalar;
44 internal::evaluator<MatrixType> matrixEval(matrix);
45 result.coeffRef(0, 0) = Scalar(1) / matrixEval.coeff(0, 0);
46 }
47};
48
49template <typename MatrixType, typename ResultType>
50struct compute_inverse_and_det_with_check<MatrixType, ResultType, 1> {
51 EIGEN_DEVICE_FUNC static inline void run(const MatrixType& matrix,
52 const typename MatrixType::RealScalar& absDeterminantThreshold,
53 ResultType& result, typename ResultType::Scalar& determinant,
54 bool& invertible) {
55 using std::abs;
56 determinant = matrix.coeff(0, 0);
57 invertible = abs(determinant) > absDeterminantThreshold;
58 if (invertible) result.coeffRef(0, 0) = typename ResultType::Scalar(1) / determinant;
59 }
60};
61
62/****************************
63*** Size 2 implementation ***
64****************************/
65
66template <typename MatrixType, typename ResultType>
67EIGEN_DEVICE_FUNC inline void compute_inverse_size2_helper(const MatrixType& matrix,
68 const typename ResultType::Scalar& invdet,
69 ResultType& result) {
70 typename ResultType::Scalar temp = matrix.coeff(0, 0);
71 result.coeffRef(0, 0) = matrix.coeff(1, 1) * invdet;
72 result.coeffRef(1, 0) = -matrix.coeff(1, 0) * invdet;
73 result.coeffRef(0, 1) = -matrix.coeff(0, 1) * invdet;
74 result.coeffRef(1, 1) = temp * invdet;
75}
76
77template <typename MatrixType, typename ResultType>
78struct compute_inverse<MatrixType, ResultType, 2> {
79 EIGEN_DEVICE_FUNC static inline void run(const MatrixType& matrix, ResultType& result) {
80 typedef typename ResultType::Scalar Scalar;
81 const Scalar invdet = typename MatrixType::Scalar(1) / matrix.determinant();
82 compute_inverse_size2_helper(matrix, invdet, result);
83 }
84};
85
86template <typename MatrixType, typename ResultType>
87struct compute_inverse_and_det_with_check<MatrixType, ResultType, 2> {
88 EIGEN_DEVICE_FUNC static inline void run(const MatrixType& matrix,
89 const typename MatrixType::RealScalar& absDeterminantThreshold,
90 ResultType& inverse, typename ResultType::Scalar& determinant,
91 bool& invertible) {
92 using std::abs;
93 typedef typename ResultType::Scalar Scalar;
94 determinant = matrix.determinant();
95 invertible = abs(determinant) > absDeterminantThreshold;
96 if (!invertible) return;
97 const Scalar invdet = Scalar(1) / determinant;
98 compute_inverse_size2_helper(matrix, invdet, inverse);
99 }
100};
101
102/****************************
103*** Size 3 implementation ***
104****************************/
105
106template <typename MatrixType, int i, int j>
107EIGEN_DEVICE_FUNC inline typename MatrixType::Scalar cofactor_3x3(const MatrixType& m) {
108 enum { i1 = (i + 1) % 3, i2 = (i + 2) % 3, j1 = (j + 1) % 3, j2 = (j + 2) % 3 };
109 return m.coeff(i1, j1) * m.coeff(i2, j2) - m.coeff(i1, j2) * m.coeff(i2, j1);
110}
111
112template <typename MatrixType, typename ResultType>
113EIGEN_DEVICE_FUNC inline void compute_inverse_size3_helper(
114 const MatrixType& matrix, const typename ResultType::Scalar& invdet,
115 const Matrix<typename ResultType::Scalar, 3, 1>& cofactors_col0, ResultType& result) {
116 // Compute cofactors in a way that avoids aliasing issues.
117 typedef typename ResultType::Scalar Scalar;
118 const Scalar c01 = cofactor_3x3<MatrixType, 0, 1>(matrix) * invdet;
119 const Scalar c11 = cofactor_3x3<MatrixType, 1, 1>(matrix) * invdet;
120 const Scalar c02 = cofactor_3x3<MatrixType, 0, 2>(matrix) * invdet;
121 result.coeffRef(1, 2) = cofactor_3x3<MatrixType, 2, 1>(matrix) * invdet;
122 result.coeffRef(2, 1) = cofactor_3x3<MatrixType, 1, 2>(matrix) * invdet;
123 result.coeffRef(2, 2) = cofactor_3x3<MatrixType, 2, 2>(matrix) * invdet;
124 result.coeffRef(1, 0) = c01;
125 result.coeffRef(1, 1) = c11;
126 result.coeffRef(2, 0) = c02;
127 result.row(0) = cofactors_col0 * invdet;
128}
129
130template <typename MatrixType, typename ResultType>
131struct compute_inverse<MatrixType, ResultType, 3> {
132 EIGEN_DEVICE_FUNC static inline void run(const MatrixType& matrix, ResultType& result) {
133 typedef typename ResultType::Scalar Scalar;
134 Matrix<typename MatrixType::Scalar, 3, 1> cofactors_col0;
135 cofactors_col0.coeffRef(0) = cofactor_3x3<MatrixType, 0, 0>(matrix);
136 cofactors_col0.coeffRef(1) = cofactor_3x3<MatrixType, 1, 0>(matrix);
137 cofactors_col0.coeffRef(2) = cofactor_3x3<MatrixType, 2, 0>(matrix);
138 const Scalar det = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
139 const Scalar invdet = Scalar(1) / det;
140 compute_inverse_size3_helper(matrix, invdet, cofactors_col0, result);
141 }
142};
143
144template <typename MatrixType, typename ResultType>
145struct compute_inverse_and_det_with_check<MatrixType, ResultType, 3> {
146 EIGEN_DEVICE_FUNC static inline void run(const MatrixType& matrix,
147 const typename MatrixType::RealScalar& absDeterminantThreshold,
148 ResultType& inverse, typename ResultType::Scalar& determinant,
149 bool& invertible) {
150 typedef typename ResultType::Scalar Scalar;
151 Matrix<Scalar, 3, 1> cofactors_col0;
152 cofactors_col0.coeffRef(0) = cofactor_3x3<MatrixType, 0, 0>(matrix);
153 cofactors_col0.coeffRef(1) = cofactor_3x3<MatrixType, 1, 0>(matrix);
154 cofactors_col0.coeffRef(2) = cofactor_3x3<MatrixType, 2, 0>(matrix);
155 determinant = (cofactors_col0.cwiseProduct(matrix.col(0))).sum();
156 invertible = Eigen::numext::abs(determinant) > absDeterminantThreshold;
157 if (!invertible) return;
158 const Scalar invdet = Scalar(1) / determinant;
159 compute_inverse_size3_helper(matrix, invdet, cofactors_col0, inverse);
160 }
161};
162
163/****************************
164*** Size 4 implementation ***
165****************************/
166
167template <typename Derived>
168EIGEN_DEVICE_FUNC inline const typename Derived::Scalar general_det3_helper(const MatrixBase<Derived>& matrix, int i1,
169 int i2, int i3, int j1, int j2, int j3) {
170 return matrix.coeff(i1, j1) *
171 (matrix.coeff(i2, j2) * matrix.coeff(i3, j3) - matrix.coeff(i2, j3) * matrix.coeff(i3, j2));
172}
173
174template <typename MatrixType, int i, int j>
175EIGEN_DEVICE_FUNC inline typename MatrixType::Scalar cofactor_4x4(const MatrixType& matrix) {
176 enum { i1 = (i + 1) % 4, i2 = (i + 2) % 4, i3 = (i + 3) % 4, j1 = (j + 1) % 4, j2 = (j + 2) % 4, j3 = (j + 3) % 4 };
177 return general_det3_helper(matrix, i1, i2, i3, j1, j2, j3) + general_det3_helper(matrix, i2, i3, i1, j1, j2, j3) +
178 general_det3_helper(matrix, i3, i1, i2, j1, j2, j3);
179}
180
181template <int Arch, typename Scalar, typename MatrixType, typename ResultType>
182struct compute_inverse_size4 {
183 EIGEN_DEVICE_FUNC static void run(const MatrixType& matrix, ResultType& result) {
184 result.coeffRef(0, 0) = cofactor_4x4<MatrixType, 0, 0>(matrix);
185 result.coeffRef(1, 0) = -cofactor_4x4<MatrixType, 0, 1>(matrix);
186 result.coeffRef(2, 0) = cofactor_4x4<MatrixType, 0, 2>(matrix);
187 result.coeffRef(3, 0) = -cofactor_4x4<MatrixType, 0, 3>(matrix);
188 result.coeffRef(0, 2) = cofactor_4x4<MatrixType, 2, 0>(matrix);
189 result.coeffRef(1, 2) = -cofactor_4x4<MatrixType, 2, 1>(matrix);
190 result.coeffRef(2, 2) = cofactor_4x4<MatrixType, 2, 2>(matrix);
191 result.coeffRef(3, 2) = -cofactor_4x4<MatrixType, 2, 3>(matrix);
192 result.coeffRef(0, 1) = -cofactor_4x4<MatrixType, 1, 0>(matrix);
193 result.coeffRef(1, 1) = cofactor_4x4<MatrixType, 1, 1>(matrix);
194 result.coeffRef(2, 1) = -cofactor_4x4<MatrixType, 1, 2>(matrix);
195 result.coeffRef(3, 1) = cofactor_4x4<MatrixType, 1, 3>(matrix);
196 result.coeffRef(0, 3) = -cofactor_4x4<MatrixType, 3, 0>(matrix);
197 result.coeffRef(1, 3) = cofactor_4x4<MatrixType, 3, 1>(matrix);
198 result.coeffRef(2, 3) = -cofactor_4x4<MatrixType, 3, 2>(matrix);
199 result.coeffRef(3, 3) = cofactor_4x4<MatrixType, 3, 3>(matrix);
200 result /= (matrix.col(0).cwiseProduct(result.row(0).transpose())).sum();
201 }
202};
203
204template <typename MatrixType, typename ResultType>
205struct compute_inverse<MatrixType, ResultType, 4>
206 : compute_inverse_size4<Architecture::Target, typename MatrixType::Scalar, MatrixType, ResultType> {};
207
208template <typename MatrixType, typename ResultType>
209struct compute_inverse_and_det_with_check<MatrixType, ResultType, 4> {
210 EIGEN_DEVICE_FUNC static inline void run(const MatrixType& matrix,
211 const typename MatrixType::RealScalar& absDeterminantThreshold,
212 ResultType& inverse, typename ResultType::Scalar& determinant,
213 bool& invertible) {
214 using std::abs;
215 determinant = matrix.determinant();
216 invertible = abs(determinant) > absDeterminantThreshold;
217 if (invertible && extract_data(matrix) != extract_data(inverse)) {
218 compute_inverse<MatrixType, ResultType>::run(matrix, inverse);
219 } else if (invertible) {
220 MatrixType matrix_t = matrix;
221 compute_inverse<MatrixType, ResultType>::run(matrix_t, inverse);
222 }
223 }
224};
225
226/*************************
227*** MatrixBase methods ***
228*************************/
229
230} // end namespace internal
231
232namespace internal {
233
234// Specialization for "dense = dense_xpr.inverse()"
235template <typename DstXprType, typename XprType>
236struct Assignment<DstXprType, Inverse<XprType>,
237 internal::assign_op<typename DstXprType::Scalar, typename XprType::Scalar>, Dense2Dense> {
238 typedef Inverse<XprType> SrcXprType;
239 EIGEN_DEVICE_FUNC static void run(DstXprType& dst, const SrcXprType& src,
240 const internal::assign_op<typename DstXprType::Scalar, typename XprType::Scalar>&) {
241 Index dstRows = src.rows();
242 Index dstCols = src.cols();
243 if ((dst.rows() != dstRows) || (dst.cols() != dstCols)) dst.resize(dstRows, dstCols);
244
245 const int Size = plain_enum_min(XprType::ColsAtCompileTime, DstXprType::ColsAtCompileTime);
246 EIGEN_ONLY_USED_FOR_DEBUG(Size);
247 eigen_assert(((Size <= 1) || (Size > 4) || (extract_data(src.nestedExpression()) != extract_data(dst))) &&
248 "Aliasing problem detected in inverse(), you need to do inverse().eval() here.");
249
250 typedef typename internal::nested_eval<XprType, XprType::ColsAtCompileTime>::type ActualXprType;
251 typedef internal::remove_all_t<ActualXprType> ActualXprTypeCleanded;
252
253 ActualXprType actual_xpr(src.nestedExpression());
254
255 compute_inverse<ActualXprTypeCleanded, DstXprType>::run(actual_xpr, dst);
256 }
257};
258
259} // end namespace internal
260
278template <typename Derived>
279EIGEN_DEVICE_FUNC inline const Inverse<Derived> MatrixBase<Derived>::inverse() const {
280 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsInteger, THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES)
281 eigen_assert(rows() == cols());
282 return Inverse<Derived>(derived());
283}
284
305template <typename Derived>
306template <typename ResultType>
307inline void MatrixBase<Derived>::computeInverseAndDetWithCheck(ResultType& inverse,
308 typename ResultType::Scalar& determinant,
309 bool& invertible,
310 const RealScalar& absDeterminantThreshold) const {
311 // i'd love to put some static assertions there, but SFINAE means that they have no effect...
312 eigen_assert(rows() == cols());
313 // for 2x2, it's worth giving a chance to avoid evaluating.
314 // for larger sizes, evaluating has negligible cost and limits code size.
315 typedef std::conditional_t<RowsAtCompileTime == 2,
316 internal::remove_all_t<typename internal::nested_eval<Derived, 2>::type>, PlainObject>
317 MatrixType;
318 internal::compute_inverse_and_det_with_check<MatrixType, ResultType>::run(derived(), absDeterminantThreshold, inverse,
319 determinant, invertible);
320}
321
341template <typename Derived>
342template <typename ResultType>
343inline void MatrixBase<Derived>::computeInverseWithCheck(ResultType& inverse, bool& invertible,
344 const RealScalar& absDeterminantThreshold) const {
345 Scalar determinant;
346 // i'd love to put some static assertions there, but SFINAE means that they have no effect...
347 eigen_assert(rows() == cols());
348 computeInverseAndDetWithCheck(inverse, determinant, invertible, absDeterminantThreshold);
349}
350
351} // end namespace Eigen
352
353#endif // EIGEN_INVERSE_IMPL_H
Eigen::DenseBase< Homogeneous< MatrixType, Direction_ > >::Scalar
internal::traits< Homogeneous< MatrixType, Direction_ > >::Scalar Scalar
Definition DenseBase.h:62
Eigen::Inverse
Expression of the inverse of another expression.
Definition Inverse.h:43
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:52
Eigen
Namespace containing all symbols from the Eigen library.
Definition Core:137
Eigen::inverse
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_inverse_op< typename Derived::Scalar >, const Derived > inverse(const Eigen::ArrayBase< Derived > &x)
Eigen::abs
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition Meta.h:523