11 #ifndef EIGEN_INVERSE_IMPL_H
12 #define EIGEN_INVERSE_IMPL_H
14 #include "./InternalHeaderCheck.h"
24 template<
typename MatrixType,
typename ResultType,
int Size = MatrixType::RowsAtCompileTime>
25 struct compute_inverse
28 static inline void run(
const MatrixType& matrix, ResultType& result)
30 result = matrix.partialPivLu().inverse();
34 template<
typename MatrixType,
typename ResultType,
int Size = MatrixType::RowsAtCompileTime>
35 struct compute_inverse_and_det_with_check { };
41 template<
typename MatrixType,
typename ResultType>
42 struct compute_inverse<MatrixType, ResultType, 1>
45 static inline void run(
const MatrixType& matrix, ResultType& result)
47 typedef typename MatrixType::Scalar Scalar;
48 internal::evaluator<MatrixType> matrixEval(matrix);
49 result.coeffRef(0,0) = Scalar(1) / matrixEval.coeff(0,0);
53 template<
typename MatrixType,
typename ResultType>
54 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 1>
57 static inline void run(
58 const MatrixType& matrix,
59 const typename MatrixType::RealScalar& absDeterminantThreshold,
61 typename ResultType::Scalar& determinant,
66 determinant = matrix.coeff(0,0);
67 invertible =
abs(determinant) > absDeterminantThreshold;
68 if(invertible) result.coeffRef(0,0) =
typename ResultType::Scalar(1) / determinant;
76 template<
typename MatrixType,
typename ResultType>
78 inline void compute_inverse_size2_helper(
79 const MatrixType& matrix,
const typename ResultType::Scalar& invdet,
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;
89 template<
typename MatrixType,
typename ResultType>
90 struct compute_inverse<MatrixType, ResultType, 2>
93 static inline void run(
const MatrixType& matrix, ResultType& result)
95 typedef typename ResultType::Scalar Scalar;
96 const Scalar invdet =
typename MatrixType::Scalar(1) / matrix.determinant();
97 compute_inverse_size2_helper(matrix, invdet, result);
101 template<
typename MatrixType,
typename ResultType>
102 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 2>
105 static inline void run(
106 const MatrixType& matrix,
107 const typename MatrixType::RealScalar& absDeterminantThreshold,
109 typename ResultType::Scalar& determinant,
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);
127 template<
typename MatrixType,
int i,
int j>
129 inline typename MatrixType::Scalar cofactor_3x3(
const MatrixType& m)
137 return m.coeff(i1, j1) * m.coeff(i2, j2)
138 - m.coeff(i1, j2) * m.coeff(i2, j1);
141 template<
typename MatrixType,
typename ResultType>
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,
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;
163 template<
typename MatrixType,
typename ResultType>
164 struct compute_inverse<MatrixType, ResultType, 3>
167 static inline void run(
const MatrixType& matrix, ResultType& result)
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);
180 template<
typename MatrixType,
typename ResultType>
181 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 3>
184 static inline void run(
185 const MatrixType& matrix,
186 const typename MatrixType::RealScalar& absDeterminantThreshold,
188 typename ResultType::Scalar& determinant,
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);
209 template<
typename Derived>
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)
214 return matrix.coeff(i1,j1)
215 * (matrix.coeff(i2,j2) * matrix.coeff(i3,j3) - matrix.coeff(i2,j3) * matrix.coeff(i3,j2));
218 template<
typename MatrixType,
int i,
int j>
220 inline typename MatrixType::Scalar cofactor_4x4(
const MatrixType& matrix)
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);
235 template<
int Arch,
typename Scalar,
typename MatrixType,
typename ResultType>
236 struct compute_inverse_size4
239 static void run(
const MatrixType& matrix, ResultType& result)
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();
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>
268 template<
typename MatrixType,
typename ResultType>
269 struct compute_inverse_and_det_with_check<MatrixType, ResultType, 4>
272 static inline void run(
273 const MatrixType& matrix,
274 const typename MatrixType::RealScalar& absDeterminantThreshold,
276 typename ResultType::Scalar& determinant,
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);
286 else if(invertible) {
287 MatrixType matrix_t = matrix;
288 compute_inverse<MatrixType, ResultType>::run(matrix_t,
inverse);
302 template<
typename DstXprType,
typename XprType>
303 struct Assignment<DstXprType, Inverse<XprType>, internal::assign_op<typename DstXprType::Scalar,typename XprType::Scalar>, Dense2Dense>
305 typedef Inverse<XprType> SrcXprType;
307 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<typename DstXprType::Scalar,typename XprType::Scalar> &)
309 Index dstRows = src.rows();
310 Index dstCols = src.cols();
311 if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
312 dst.resize(dstRows, dstCols);
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.");
319 typedef typename internal::nested_eval<XprType,XprType::ColsAtCompileTime>::type ActualXprType;
320 typedef internal::remove_all_t<ActualXprType> ActualXprTypeCleanded;
322 ActualXprType actual_xpr(src.nestedExpression());
324 compute_inverse<ActualXprTypeCleanded, DstXprType>::run(actual_xpr, dst);
348 template<
typename Derived>
353 eigen_assert(rows() == cols());
377 template<
typename Derived>
378 template<
typename ResultType>
381 typename ResultType::Scalar& determinant,
383 const RealScalar& absDeterminantThreshold
387 eigen_assert(rows() == cols());
390 typedef std::conditional_t<
391 RowsAtCompileTime == 2,
392 internal::remove_all_t<typename internal::nested_eval<Derived, 2>::type>,
395 internal::compute_inverse_and_det_with_check<MatrixType, ResultType>::run
396 (derived(), absDeterminantThreshold,
inverse, determinant, invertible);
418 template<
typename Derived>
419 template<
typename ResultType>
423 const RealScalar& absDeterminantThreshold
428 eigen_assert(rows() == cols());
429 computeInverseAndDetWithCheck(
inverse,determinant,invertible,absDeterminantThreshold);
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