1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_INVERSE_H
#define EIGEN_INVERSE_H
/********************************************************************
*** Part 1 : optimized implementations for fixed-size 2,3,4 cases ***
********************************************************************/
template<typename MatrixType>
void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar invdet = Scalar(1) / matrix.determinant();
result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
result->coeffRef(1,1) = matrix.coeff(0,0) * invdet;
}
template<typename XprType, typename MatrixType>
bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar det = matrix.determinant();
if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false;
const Scalar invdet = Scalar(1) / det;
result->coeffRef(0,0) = matrix.coeff(1,1) * invdet;
result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
result->coeffRef(1,1) = matrix.coeff(0,0) * invdet;
return true;
}
template<typename MatrixType>
void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
const Scalar det_minor00 = matrix.minor(0,0).determinant();
const Scalar det_minor10 = matrix.minor(1,0).determinant();
const Scalar det_minor20 = matrix.minor(2,0).determinant();
const Scalar invdet = Scalar(1) / ( det_minor00 * matrix.coeff(0,0)
- det_minor10 * matrix.coeff(1,0)
+ det_minor20 * matrix.coeff(2,0) );
result->coeffRef(0, 0) = det_minor00 * invdet;
result->coeffRef(0, 1) = -det_minor10 * invdet;
result->coeffRef(0, 2) = det_minor20 * invdet;
result->coeffRef(1, 0) = -matrix.minor(0,1).determinant() * invdet;
result->coeffRef(1, 1) = matrix.minor(1,1).determinant() * invdet;
result->coeffRef(1, 2) = -matrix.minor(2,1).determinant() * invdet;
result->coeffRef(2, 0) = matrix.minor(0,2).determinant() * invdet;
result->coeffRef(2, 1) = -matrix.minor(1,2).determinant() * invdet;
result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
}
template<typename MatrixType>
bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixType* result)
{
/* Let's split M into four 2x2 blocks:
* (P Q)
* (R S)
* If P is invertible, with inverse denoted by P_inverse, and if
* (S - R*P_inverse*Q) is also invertible, then the inverse of M is
* (P' Q')
* (R' S')
* where
* S' = (S - R*P_inverse*Q)^(-1)
* P' = P1 + (P1*Q) * S' *(R*P_inverse)
* Q' = -(P_inverse*Q) * S'
* R' = -S' * (R*P_inverse)
*/
typedef Block<MatrixType,2,2> XprBlock22;
typedef typename XprBlock22::Eval Block22;
Block22 P_inverse;
if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse))
{
const Block22 Q = matrix.template block<2,2>(0,2);
const Block22 P_inverse_times_Q = P_inverse * Q;
const XprBlock22 R = matrix.template block<2,2>(2,0);
const Block22 R_times_P_inverse = R * P_inverse;
const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q;
const XprBlock22 S = matrix.template block<2,2>(2,2);
const Block22 X = S - R_times_P_inverse_times_Q;
Block22 Y;
ei_compute_inverse_in_size2_case(X, &Y);
result->template block<2,2>(2,2) = Y;
result->template block<2,2>(2,0) = - Y * R_times_P_inverse;
const Block22 Z = P_inverse_times_Q * Y;
result->template block<2,2>(0,2) = - Z;
result->template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse;
return true;
}
else
{
return false;
}
}
template<typename MatrixType>
void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result)
{
if(ei_compute_inverse_in_size4_case_helper(matrix, result))
{
// good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful.
return;
}
else
{
// rare case: the topleft 2x2 block is not invertible (but the matrix itself is assumed to be).
// since this is a rare case, we don't need to optimize it. We just want to handle it with little
// additional code.
MatrixType m(matrix);
m.row(1).swap(m.row(2));
if(ei_compute_inverse_in_size4_case_helper(m, result))
{
// good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that two
// rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting
// the corresponding columns.
result->col(1).swap(result->col(2));
}
else
{
// last possible case. Since matrix is assumed to be invertible, this last case has to work.
m.row(1).swap(m.row(2));
m.row(1).swap(m.row(3));
ei_compute_inverse_in_size4_case_helper(m, result);
result->col(1).swap(result->col(3));
}
}
}
/***********************************************
*** Part 2 : selector and MatrixBase methods ***
***********************************************/
template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
struct ei_compute_inverse
{
static inline void run(const MatrixType& matrix, MatrixType* result)
{
LU<MatrixType> lu(matrix);
lu.computeInverse(result);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 1>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
{
typedef typename MatrixType::Scalar Scalar;
result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 2>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
{
ei_compute_inverse_in_size2_case(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 3>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
{
ei_compute_inverse_in_size3_case(matrix, result);
}
};
template<typename MatrixType>
struct ei_compute_inverse<MatrixType, 4>
{
static inline void run(const MatrixType& matrix, MatrixType* result)
{
ei_compute_inverse_in_size4_case(matrix, result);
}
};
/** \lu_module
*
* Computes the matrix inverse of this matrix.
*
* \note This matrix must be invertible, otherwise the result is undefined.
*
* \param result Pointer to the matrix in which to store the result.
*
* Example: \include MatrixBase_computeInverse.cpp
* Output: \verbinclude MatrixBase_computeInverse.out
*
* \sa inverse()
*/
template<typename Derived>
inline void MatrixBase<Derived>::computeInverse(EvalType *result) const
{
typedef typename ei_eval<Derived>::type MatrixType;
ei_assert(rows() == cols());
EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,scalar_type_must_be_floating_point);
ei_compute_inverse<MatrixType>::run(eval(), result);
}
/** \lu_module
*
* \returns the matrix inverse of this matrix.
*
* \note This matrix must be invertible, otherwise the result is undefined.
*
* \note This method returns a matrix by value, which can be inefficient. To avoid that overhead,
* use computeInverse() instead.
*
* Example: \include MatrixBase_inverse.cpp
* Output: \verbinclude MatrixBase_inverse.out
*
* \sa computeInverse()
*/
template<typename Derived>
inline const typename MatrixBase<Derived>::EvalType MatrixBase<Derived>::inverse() const
{
EvalType result(rows(), cols());
computeInverse(&result);
return result;
}
#endif // EIGEN_INVERSE_H
|
