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