STK++ 0.9.13
STK_Invert.h
Go to the documentation of this file.
1/*--------------------------------------------------------------------*/
2/* Copyright (C) 2004-2016 Serge Iovleff, Université Lille 1, Inria
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU Lesser General Public License as
6 published by the Free Software Foundation; either version 2 of the
7 License, or (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU Lesser General Public License for more details.
13
14 You should have received a copy of the GNU Lesser General Public
15 License along with this program; if not, write to the
16 Free Software Foundation, Inc.,
17 59 Temple Place,
18 Suite 330,
19 Boston, MA 02111-1307
20 USA
21
22 Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
23*/
24
25/*
26 * Project: stkpp::Algebra
27 * created on: 9 août 2016
28 * Author: iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
29 **/
30
35#ifndef STK_INVERT_H
36#define STK_INVERT_H
37
38#ifdef STKUSELAPACK
39#include <Algebra/include/STK_lapack_Svd.h>
40#else
41#include <Algebra/include/STK_Svd.h>
42#endif
43
44namespace STK
45{
46
47namespace hidden
48{
52template<class Matrix, int Size_>
53struct InvertMatrixImpl
54{
55 typedef typename Matrix::Type Type;
61 static Type invertMatrix11( Matrix const& m, CArraySquare<Type, Size_>& inv)
62 {
63 const int iShift = m.beginRows(), jShift = m.beginCols();
64 inv.resize(TRange<Size_>(0, 1));
65 // cofactor (0, jShift+0) [0]
66 inv(0, 0) = Type(1);
67 // compute determinant
68 Type det = m(iShift+0, jShift+0);
69 if (det == Type(0)) return Type(0);
70 // compute inverse matrix
71 inv /= det;
72 return det;
73 }
79 static Type invertMatrix22( Matrix const& m, CArraySquare<Type, Size_>& inv)
80 {
81 const int iShift = m.beginRows(), jShift = m.beginCols();
82 inv.resize(TRange<Size_>(0, 2));
83 // cofactor (0, jShift+0) [0]
84 inv(0, 0) = m(iShift+1, jShift+1);
85 // cofactor (1, jShift+0) [1]
86 inv(1, 0) = - m(iShift+1, jShift+0);
87 // cofactor (1, jShift+1) [2]
88 inv(0, 1) = - m(iShift+0, jShift+1);
89 // cofactor (1, jShift+1) [3]
90 inv(1, 1) = m(iShift+0, jShift+0);
91
92 // determinant
93 Type det = m(iShift+0, jShift+0) * inv(0, 0)
94 + m(iShift+0, jShift+1) * inv(1, 0);
95 if (det == Type(0)) return Type(0);
96 // compute inverse matrix
97 inv /= det;
98 return det;
99 }
105 static Type invertMatrix33( Matrix const& m, CArraySquare<Type, Size_>& inv)
106 {
107 const int iShift = m.beginRows(), jShift = m.beginCols();
108 inv.resize(TRange<Size_>(0, 3));
109 // cofactor
110 inv(0, 0) = m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
111 - m(iShift+2, jShift+1) * m(iShift+1, jShift+2);
112 inv(0, 1) = m(iShift+0, jShift+2) * m(iShift+2, jShift+1)
113 - m(iShift+0, jShift+1) * m(iShift+2, jShift+2);
114 inv(0, 2) = m(iShift+0, jShift+1) * m(iShift+1, jShift+2)
115 - m(iShift+0, jShift+2) * m(iShift+1, jShift+1);
116 inv(1, 0) = m(iShift+1, jShift+2) * m(iShift+2, jShift+0)
117 - m(iShift+1, jShift+0) * m(iShift+2, jShift+2);
118 inv(1, 1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+2)
119 - m(iShift+0, jShift+2) * m(iShift+2, jShift+0);
120 inv(1, 2) = m(iShift+1, jShift+0) * m(iShift+0, jShift+2)
121 - m(iShift+0, jShift+0) * m(iShift+1, jShift+2);
122 inv(2, 0) = m(iShift+1, jShift+0) * m(iShift+2, jShift+1)
123 - m(iShift+2, jShift+0) * m(iShift+1, jShift+1);
124 inv(2, 1) = m(iShift+2, jShift+0) * m(iShift+0, jShift+1)
125 - m(iShift+0, jShift+0) * m(iShift+2, jShift+1);
126 inv(2, 2) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1)
127 - m(iShift+1, jShift+0) * m(iShift+0, jShift+1);
128
129 // computes the inverse of a matrix m
130 Type det = m(iShift+0, jShift+0) * inv(0, 0)
131 + m(iShift+0, jShift+1) * inv(1, 0)
132 + m(iShift+0, jShift+2) * inv(2, 0);
133 if (det == Type(0)) return Type(0);
134 // compute inverse matrix
135 inv /= det;
136 return det;
137 }
143 static Type invertMatrix44( Matrix const& m, CArraySquare<Type, Size_>& inv)
144 {
145 const int iShift = m.beginRows(), jShift = m.beginCols();
146 inv.resize(TRange<Size_>(0, 4));
147 // cofactor (0, 0) [0]
148 inv(0, 0) = m(iShift+1, jShift+1) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
149 - m(iShift+1, jShift+1) * m(iShift+3, jShift+2) * m(iShift+2, jShift+3)
150 - m(iShift+1, jShift+2) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
151 + m(iShift+1, jShift+2) * m(iShift+3, jShift+1) * m(iShift+2, jShift+3)
152 + m(iShift+1, jShift+3) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
153 - m(iShift+1, jShift+3) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2);
154 // cofactor (0, 1) [4]
155 inv(0, 1) = -m(iShift+0, jShift+1) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
156 + m(iShift+0, jShift+1) * m(iShift+3, jShift+2) * m(iShift+2, jShift+3)
157 + m(iShift+0, jShift+2) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
158 - m(iShift+0, jShift+2) * m(iShift+3, jShift+1) * m(iShift+2, jShift+3)
159 - m(iShift+0, jShift+3) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
160 + m(iShift+0, jShift+3) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2);
161 // cofactor (0, 2) [8]
162 inv(0, 2) = m(iShift+0, jShift+1) * m(iShift+1, jShift+2) * m(iShift+3, jShift+3)
163 - m(iShift+0, jShift+1) * m(iShift+3, jShift+2) * m(iShift+1, jShift+3)
164 - m(iShift+0, jShift+2) * m(iShift+1, jShift+1) * m(iShift+3, jShift+3)
165 + m(iShift+0, jShift+2) * m(iShift+3, jShift+1) * m(iShift+1, jShift+3)
166 + m(iShift+0, jShift+3) * m(iShift+1, jShift+1) * m(iShift+3, jShift+2)
167 - m(iShift+0, jShift+3) * m(iShift+3, jShift+1) * m(iShift+1, jShift+2);
168 // cofactor (0, 3) [12]
169 inv(0, 3) = -m(iShift+0, jShift+1) * m(iShift+1, jShift+2) * m(iShift+2, jShift+3)
170 + m(iShift+0, jShift+1) * m(iShift+2, jShift+2) * m(iShift+1, jShift+3)
171 + m(iShift+0, jShift+2) * m(iShift+1, jShift+1) * m(iShift+2, jShift+3)
172 - m(iShift+0, jShift+2) * m(iShift+2, jShift+1) * m(iShift+1, jShift+3)
173 - m(iShift+0, jShift+3) * m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
174 + m(iShift+0, jShift+3) * m(iShift+2, jShift+1) * m(iShift+1, jShift+2);
175 // cofactor (1, 0) [1]
176 inv(1, 0) = -m(iShift+1, jShift+0) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
177 + m(iShift+1, jShift+0) * m(iShift+3, jShift+2) * m(iShift+2, jShift+3)
178 + m(iShift+1, jShift+2) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
179 - m(iShift+1, jShift+2) * m(iShift+3, jShift+0) * m(iShift+2, jShift+3)
180 - m(iShift+1, jShift+3) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
181 + m(iShift+1, jShift+3) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2);
182 // cofactor (1, 1) [5]
183 inv(1, 1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
184 - m(iShift+0, jShift+0) * m(iShift+3, jShift+2) * m(iShift+2, jShift+3)
185 - m(iShift+0, jShift+2) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
186 + m(iShift+0, jShift+2) * m(iShift+3, jShift+0) * m(iShift+2, jShift+3)
187 + m(iShift+0, jShift+3) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
188 - m(iShift+0, jShift+3) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2);
189 // cofactor (1, 2) [9]
190 inv(1, 2) = -m(iShift+0, jShift+0) * m(iShift+1, jShift+2) * m(iShift+3, jShift+3)
191 + m(iShift+0, jShift+0) * m(iShift+3, jShift+2) * m(iShift+1, jShift+3)
192 + m(iShift+0, jShift+2) * m(iShift+1, jShift+0) * m(iShift+3, jShift+3)
193 - m(iShift+0, jShift+2) * m(iShift+3, jShift+0) * m(iShift+1, jShift+3)
194 - m(iShift+0, jShift+3) * m(iShift+1, jShift+0) * m(iShift+3, jShift+2)
195 + m(iShift+0, jShift+3) * m(iShift+3, jShift+0) * m(iShift+1, jShift+2);
196 // cofactor (1, 3) [13]
197 inv(1, 3) = m(iShift+0, jShift+0) * m(iShift+1, jShift+2) * m(iShift+2, jShift+3)
198 - m(iShift+0, jShift+0) * m(iShift+2, jShift+2) * m(iShift+1, jShift+3)
199 - m(iShift+0, jShift+2) * m(iShift+1, jShift+0) * m(iShift+2, jShift+3)
200 + m(iShift+0, jShift+2) * m(iShift+2, jShift+0) * m(iShift+1, jShift+3)
201 + m(iShift+0, jShift+3) * m(iShift+1, jShift+0) * m(iShift+2, jShift+2)
202 - m(iShift+0, jShift+3) * m(iShift+2, jShift+0) * m(iShift+1, jShift+2);
203 // cofactor (2, 0) [2]
204 inv(2, 0) = m(iShift+1, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
205 - m(iShift+1, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+3)
206 - m(iShift+1, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
207 + m(iShift+1, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+3)
208 + m(iShift+1, jShift+3) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
209 - m(iShift+1, jShift+3) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
210 // cofactor (2, 1)
211 inv(2, 1) = -m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
212 + m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+3)
213 + m(iShift+0, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
214 - m(iShift+0, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+3)
215 - m(iShift+0, jShift+3) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
216 + m(iShift+0, jShift+3) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
217 // cofactor (2, 2)
218 inv(2, 2) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+3, jShift+3)
219 - m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+1, jShift+3)
220 - m(iShift+0, jShift+1) * m(iShift+1, jShift+0) * m(iShift+3, jShift+3)
221 + m(iShift+0, jShift+1) * m(iShift+3, jShift+0) * m(iShift+1, jShift+3)
222 + m(iShift+0, jShift+3) * m(iShift+1, jShift+0) * m(iShift+3, jShift+1)
223 - m(iShift+0, jShift+3) * m(iShift+3, jShift+0) * m(iShift+1, jShift+1);
224 // cofactor (2, 3)
225 inv(2, 3) = -m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+2, jShift+3)
226 + m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+1, jShift+3)
227 + m(iShift+0, jShift+1) * m(iShift+1, jShift+0) * m(iShift+2, jShift+3)
228 - m(iShift+0, jShift+1) * m(iShift+2, jShift+0) * m(iShift+1, jShift+3)
229 - m(iShift+0, jShift+3) * m(iShift+1, jShift+0) * m(iShift+2, jShift+1)
230 + m(iShift+0, jShift+3) * m(iShift+2, jShift+0) * m(iShift+1, jShift+1);
231 // cofactor (3, 0)
232 inv(3, 0) = -m(iShift+1, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
233 + m(iShift+1, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2)
234 + m(iShift+1, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
235 - m(iShift+1, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2)
236 - m(iShift+1, jShift+2) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
237 + m(iShift+1, jShift+2) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
238 // cofactor (3, 1) [7]
239 inv(3, 1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
240 - m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2)
241 - m(iShift+0, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
242 + m(iShift+0, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2)
243 + m(iShift+0, jShift+2) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
244 - m(iShift+0, jShift+2) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
245 // cofactor (3, 2) [11]
246 inv(3, 2) = -m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+3, jShift+2)
247 + m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+1, jShift+2)
248 + m(iShift+0, jShift+1) * m(iShift+1, jShift+0) * m(iShift+3, jShift+2)
249 - m(iShift+0, jShift+1) * m(iShift+3, jShift+0) * m(iShift+1, jShift+2)
250 - m(iShift+0, jShift+2) * m(iShift+1, jShift+0) * m(iShift+3, jShift+1)
251 + m(iShift+0, jShift+2) * m(iShift+3, jShift+0) * m(iShift+1, jShift+1);
252 // cofactor (3, 3) [15]
253 inv(3, 3) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
254 - m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+1, jShift+2)
255 - m(iShift+0, jShift+1) * m(iShift+1, jShift+0) * m(iShift+2, jShift+2)
256 + m(iShift+0, jShift+1) * m(iShift+2, jShift+0) * m(iShift+1, jShift+2)
257 + m(iShift+0, jShift+2) * m(iShift+1, jShift+0) * m(iShift+2, jShift+1)
258 - m(iShift+0, jShift+2) * m(iShift+2, jShift+0) * m(iShift+1, jShift+1);
259 // compute determinant and inverse matrix
260 Type det = m(iShift+0, jShift+0) * inv(0, 0)
261 + m(iShift+1, jShift+0) * inv(0, 1)
262 + m(iShift+2, jShift+0) * inv(0, 2)
263 + m(iShift+3, jShift+0) * inv(0, 3);
264 if (det == Type(0)) return Type(0);
265 // compute inverse matrix
266 inv /= det;
267 return det;
268 }
269
276 static Type invertMatrixXX( Matrix const& m, CArraySquare<Type, Size_>& inv)
277 {
278 inv.resize(TRange<Size_>(0, m.sizeRows()));
279#ifdef STKUSELAPACK
280 lapack::Svd decomp(m);
281#else
282 Svd<Matrix> decomp(m);
283#endif
284 if (!decomp.run()) return Type(0);
285 // compute (generalized) inverse matrix
286 decomp.ginv(inv);
287 return decomp.det();
288 }
289};
290
291
294template<class Matrix, int Size_>
295struct InvertMatrixDispatcher
296{
297 typedef typename Matrix::Type Type;
298 inline static Type run( Matrix const& m, CArraySquare<Type, Size_>& inv)
299 {
300 switch (m.sizeRows())
301 {
302 case 1: return InvertMatrixImpl<Matrix, Size_>::invertMatrix11(m, inv);
303 case 2: return InvertMatrixImpl<Matrix, Size_>::invertMatrix22(m, inv);
304 case 3: return InvertMatrixImpl<Matrix, Size_>::invertMatrix33(m, inv);
305 case 4: return InvertMatrixImpl<Matrix, Size_>::invertMatrix44(m, inv);
306 default:return InvertMatrixImpl<Matrix, Size_>::invertMatrixXX(m, inv);
307 }
308 }
309};
310
311
312template<class Matrix>
313struct InvertMatrixDispatcher<Matrix, 1>
314{
315 typedef typename Matrix::Type Type;
316 inline static Type run( Matrix const& m, CArraySquare<Type, 1>& inv)
317 { return InvertMatrixImpl<Matrix, 1>::invertMatrix11(m, inv);}
318};
319
320template<class Matrix>
321struct InvertMatrixDispatcher<Matrix, 2>
322{
323 typedef typename Matrix::Type Type;
324 inline static Type run( Matrix const& m, CArraySquare<Type, 2>& inv)
325 { return InvertMatrixImpl<Matrix, 2>::invertMatrix22(m, inv);}
326};
327
328template<class Matrix>
329struct InvertMatrixDispatcher<Matrix, 3>
330{
331 typedef typename Matrix::Type Type;
332 inline static Type run( Matrix const& m, CArraySquare<Type, 3>& inv)
333 { return InvertMatrixImpl<Matrix, 3>::invertMatrix33(m, inv);}
334};
335
336template<class Matrix>
337struct InvertMatrixDispatcher<Matrix, 4>
338{
339 typedef typename Matrix::Type Type;
340 inline static Type run( Matrix const& m, CArraySquare<Type, 4>& inv)
341 { return InvertMatrixImpl<Matrix, 4>::invertMatrix44(m, inv);}
342};
343
344} // namespace hidden
345
346} // namespace STK
347
348#endif /* STK_INVERT_H */
STK_Svd.h
In this file we define the Svd Class.
STK_lapack_Svd.h
In this file we define the enclosing class of the dgeqrf lapack routine.
STK::ISvd::run
virtual bool run()
implement the run method
Definition STK_ISvd.h:155
STK::ISvd::det
Type det() const
Definition STK_ISvd.h:140
STK::ISvd::ginv
OtherArray & ginv(OtherArray &res) const
Compute the generalized inverse of the matrix and put the result in res.
Definition STK_ISvd.h:227
STK::MultidimRegression
The MultidimRegression class allows to regress a multidimensional output variable among a multivariat...
Definition STK_MultidimRegression.h:52
STK::hidden::InvertMatrixImpl::invertMatrix44
static Type invertMatrix44(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the 4x4 matrix m and store the result in inv.
Definition STK_Invert.h:143
STK::hidden::InvertMatrixImpl::invertMatrix33
static Type invertMatrix33(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the 3x3 matrix m and store the result in inv.
Definition STK_Invert.h:105
STK::hidden::InvertMatrixImpl::invertMatrix11
static Type invertMatrix11(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the 1x1 matrix and store the result in inv.
Definition STK_Invert.h:61
STK::hidden::InvertMatrixImpl::invertMatrix22
static Type invertMatrix22(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the 2x2 matrix m and store the result in inv.
Definition STK_Invert.h:79
STK::hidden::InvertMatrixImpl::invertMatrixXX
static Type invertMatrixXX(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the matrix m and store the result in inv.
Definition STK_Invert.h:276
STK
The namespace STK is the main domain space of the Statistical ToolKit project.
STK::hidden::InvertMatrixDispatcher< Matrix, 1 >::run
static Type run(Matrix const &m, CArraySquare< Type, 1 > &inv)
Definition STK_Invert.h:316
STK::hidden::InvertMatrixDispatcher< Matrix, 2 >::run
static Type run(Matrix const &m, CArraySquare< Type, 2 > &inv)
Definition STK_Invert.h:324
STK::hidden::InvertMatrixDispatcher< Matrix, 3 >::run
static Type run(Matrix const &m, CArraySquare< Type, 3 > &inv)
Definition STK_Invert.h:332
STK::hidden::InvertMatrixDispatcher< Matrix, 4 >::run
static Type run(Matrix const &m, CArraySquare< Type, 4 > &inv)
Definition STK_Invert.h:340
STK::hidden::InvertMatrixDispatcher
computing the inverse of a matrix.
Definition STK_Invert.h:296
STK::hidden::InvertMatrixDispatcher::run
static Type run(Matrix const &m, CArraySquare< Type, Size_ > &inv)
Definition STK_Invert.h:298
STK::hidden::InvertMatrixImpl
Implementation of the inversion matrix method for general/square matrices.
Definition STK_Invert.h:54