STK++ 0.9.13
STK_InvertUpperSym.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_UPPERSYMINVERT_H
36#define STK_UPPERSYMINVERT_H
37
38#ifdef STKUSELAPACK
39#include <Algebra/include/STK_lapack_SymEigen.h>
40#else
41#include <Algebra/include/STK_SymEigen.h>
42#endif
43
44namespace STK
45{
46
47namespace hidden
48{
52template<class Matrix, int Size_>
53struct InvertUpperSymMatrixImpl
54{
55 typedef typename Matrix::Type Type;
62 static Type invertUpperSymMatrix11( Matrix const& m, CArraySquare<Type, Size_>& inv)
63 {
64 const int iShift = m.beginRows(), jShift = m.beginCols();
65 inv.resize(TRange<Size_>(0, 1));
66 // cofactor (0,0) [0]
67 inv(0, 0) = Type(1);
68 // compute determinant
69 Type det = m(iShift+0, jShift+0);
70 if (det == Type(0)) return Type(0);
71 // compute inverse matrix
72 inv /= det;
73 return det;
74 }
82 static Type invertUpperSymMatrix22( Matrix const& m, CArraySquare<Type, Size_>& inv)
83 {
84 const int iShift = m.beginRows(), jShift = m.beginCols();
85 inv.resize(TRange<Size_>(0, 2));
86 // cofactor (0,0) [0]
87 inv(0, 0) = m(iShift+1, jShift+1);
88 // cofactor (1,0) [1]
89 inv(0, 1) = - m(iShift+0, jShift+1);
90 // cofactor (1,1) [3]
91 inv(1, 1) = m(iShift+0, jShift+0);
92 // symmetry
93 inv(1, 0) = inv(0, 1);
94 // determinant
95 Type det = m(iShift+0, jShift+0) * inv(0, 0)
96 + m(iShift+0, jShift+1) * inv(0, 1);
97 if (det == Type(0)) return Type(0);
98 // inverse matrix
99 inv /= det;
100 return det;
101 }
108 static Type invertUpperSymMatrix33( Matrix const& m, CArraySquare<Type, Size_>& inv)
109 {
110 const int iShift = m.beginRows(), jShift = m.beginCols();
111 inv.resize(TRange<Size_>(0, 3));
112 // cofactor (0,0) [0]
113 inv(0, 0) = m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
114 - m(iShift+1, jShift+2) * m(iShift+1, jShift+2);
115 // cofactor (1,0) [1]
116 inv(0, 1) = m(iShift+1, jShift+2) * m(iShift+0, jShift+2)
117 - m(iShift+0, jShift+1) * m(iShift+2, jShift+2);
118 inv(0, 2) = m(iShift+0, jShift+1) * m(iShift+1, jShift+2)
119 - m(iShift+0, jShift+2) * m(iShift+1, jShift+1);
120 inv(1, 1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+2)
121 - m(iShift+0, jShift+2) * m(iShift+0, jShift+2);
122 inv(1, 2) = m(iShift+0, jShift+2) * m(iShift+0, jShift+1)
123 - m(iShift+0, jShift+0) * m(iShift+1, jShift+2);
124 inv(2, 2) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1)
125 - m(iShift+0, jShift+1) * m(iShift+0, jShift+1);
126 // solution symmetric
127 inv(1, 0) = inv(0, 1); inv(2, 0) = inv(0, 2); inv(2, 1) = inv(1, 2);
128 // compute determinant and inverse matrix
129 Type det = m(iShift+0, jShift+0) * inv(0, 0)
130 + m(iShift+0, jShift+1) * inv(0, 1)
131 + m(iShift+0, jShift+2) * inv(0, 2);
132 if (det == Type(0)) return Type(0);
133 // inverse matrix
134 inv /= det;
135 return det;
136 }
143 static Type invertUpperSymMatrix44( 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+2, jShift+3) * m(iShift+2, jShift+3)
150 - m(iShift+1, jShift+2) * m(iShift+1, jShift+2) * m(iShift+3, jShift+3)
151 + m(iShift+1, jShift+2) * m(iShift+1, jShift+3) * m(iShift+2, jShift+3)
152 + m(iShift+1, jShift+3) * m(iShift+1, jShift+2) * m(iShift+2, jShift+3)
153 - m(iShift+1, jShift+3) * m(iShift+1, jShift+3) * m(iShift+2, jShift+2);
154 // cofactor (1,0) [1]
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+2, jShift+3) * m(iShift+2, jShift+3)
157 + m(iShift+1, jShift+2) * m(iShift+0, jShift+2) * m(iShift+3, jShift+3)
158 - m(iShift+1, jShift+2) * m(iShift+0, jShift+3) * m(iShift+2, jShift+3)
159 - m(iShift+1, jShift+3) * m(iShift+0, jShift+2) * m(iShift+2, jShift+3)
160 + m(iShift+1, jShift+3) * m(iShift+0, jShift+3) * m(iShift+2, jShift+2);
161 // cofactor (2,0) [2]
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+1, jShift+3) * m(iShift+2, jShift+3)
164 - m(iShift+1, jShift+1) * m(iShift+0, jShift+2) * m(iShift+3, jShift+3)
165 + m(iShift+1, jShift+1) * m(iShift+0, jShift+3) * m(iShift+2, jShift+3)
166 + m(iShift+1, jShift+3) * m(iShift+0, jShift+2) * m(iShift+1, jShift+3)
167 - m(iShift+1, jShift+3) * m(iShift+0, jShift+3) * m(iShift+1, jShift+2);
168 // cofactor (3,0) [3]
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+1, jShift+3) * m(iShift+2, jShift+2)
171 + m(iShift+1, jShift+1) * m(iShift+0, jShift+2) * m(iShift+2, jShift+3)
172 - m(iShift+1, jShift+1) * m(iShift+0, jShift+3) * m(iShift+2, jShift+2)
173 - m(iShift+1, jShift+2) * m(iShift+0, jShift+2) * m(iShift+1, jShift+3)
174 + m(iShift+1, jShift+2) * m(iShift+0, jShift+3) * m(iShift+1, jShift+2);
175 // cofactor (1,1) [5]
176 inv(1, 1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
177 - m(iShift+0, jShift+0) * m(iShift+2, jShift+3) * m(iShift+2, jShift+3)
178 - m(iShift+0, jShift+2) * m(iShift+0, jShift+2) * m(iShift+3, jShift+3)
179 + m(iShift+0, jShift+2) * m(iShift+0, jShift+3) * m(iShift+2, jShift+3)
180 + m(iShift+0, jShift+3) * m(iShift+0, jShift+2) * m(iShift+2, jShift+3)
181 - m(iShift+0, jShift+3) * m(iShift+0, jShift+3) * m(iShift+2, jShift+2);
182 // cofactor (2,1) [6]
183 inv(1, 2) = -m(iShift+0, jShift+0) * m(iShift+1, jShift+2) * m(iShift+3, jShift+3)
184 + m(iShift+0, jShift+0) * m(iShift+1, jShift+3) * m(iShift+2, jShift+3)
185 + m(iShift+0, jShift+1) * m(iShift+0, jShift+2) * m(iShift+3, jShift+3)
186 - m(iShift+0, jShift+1) * m(iShift+0, jShift+3) * m(iShift+2, jShift+3)
187 - m(iShift+0, jShift+3) * m(iShift+0, jShift+2) * m(iShift+1, jShift+3)
188 + m(iShift+0, jShift+3) * m(iShift+0, jShift+3) * m(iShift+1, jShift+2);
189 // cofactor (3,1) [7]
190 inv(1, 3) = m(iShift+0, jShift+0) * m(iShift+1, jShift+2) * m(iShift+2, jShift+3)
191 - m(iShift+0, jShift+0) * m(iShift+1, jShift+3) * m(iShift+2, jShift+2)
192 - m(iShift+0, jShift+1) * m(iShift+0, jShift+2) * m(iShift+2, jShift+3)
193 + m(iShift+0, jShift+1) * m(iShift+0, jShift+3) * m(iShift+2, jShift+2)
194 + m(iShift+0, jShift+2) * m(iShift+0, jShift+2) * m(iShift+1, jShift+3)
195 - m(iShift+0, jShift+2) * m(iShift+0, jShift+3) * m(iShift+1, jShift+2);
196 // cofactor (2,2) [10]
197 inv(2, 2) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+3, jShift+3)
198 - m(iShift+0, jShift+0) * m(iShift+1, jShift+3) * m(iShift+1, jShift+3)
199 - m(iShift+0, jShift+1) * m(iShift+0, jShift+1) * m(iShift+3, jShift+3)
200 + m(iShift+0, jShift+1) * m(iShift+0, jShift+3) * m(iShift+1, jShift+3)
201 + m(iShift+0, jShift+3) * m(iShift+0, jShift+1) * m(iShift+1, jShift+3)
202 - m(iShift+0, jShift+3) * m(iShift+0, jShift+3) * m(iShift+1, jShift+1);
203 // cofactor (3,2) [11]
204 inv(2, 3) = -m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+2, jShift+3)
205 + m(iShift+0, jShift+0) * m(iShift+1, jShift+3) * m(iShift+1, jShift+2)
206 + m(iShift+0, jShift+1) * m(iShift+0, jShift+1) * m(iShift+2, jShift+3)
207 - m(iShift+0, jShift+1) * m(iShift+0, jShift+3) * m(iShift+1, jShift+2)
208 - m(iShift+0, jShift+2) * m(iShift+0, jShift+1) * m(iShift+1, jShift+3)
209 + m(iShift+0, jShift+2) * m(iShift+0, jShift+3) * m(iShift+1, jShift+1);
210 // cofactor (3,3) [15]
211 inv(3, 3) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
212 - m(iShift+0, jShift+0) * m(iShift+1, jShift+2) * m(iShift+1, jShift+2)
213 - m(iShift+0, jShift+1) * m(iShift+0, jShift+1) * m(iShift+2, jShift+2)
214 + m(iShift+0, jShift+1) * m(iShift+0, jShift+2) * m(iShift+1, jShift+2)
215 + m(iShift+0, jShift+2) * m(iShift+0, jShift+1) * m(iShift+1, jShift+2)
216 - m(iShift+0, jShift+2) * m(iShift+0, jShift+2) * m(iShift+1, jShift+1);
217 // symmetric
218 inv(1, 0) = inv(0, 1); inv(2, 0) = inv(0, 2); inv(3, 0) = inv(0, 3);
219 inv(2, 1) = inv(1, 2); inv(3, 1) = inv(1, 3);
220 inv(3, 2) = inv(2, 3);
221 // compute determinant
222 Type det = m(iShift+0, jShift+0) * inv(0, 0)
223 + m(iShift+0, jShift+1) * inv(0, 1)
224 + m(iShift+0, jShift+2) * inv(0, 2)
225 + m(iShift+0, jShift+3) * inv(0, 3);
226 if (det == Type(0)) return Type(0);
227 // inverse matrix
228 inv /= det;
229 return det;
230 }
238 static Type invertUpperSymMatrixXX( Matrix const& m, CArraySquare<Type, Size_>& inv)
239 {
240#ifdef STKUSELAPACK
241 lapack::SymEigen<Matrix> decomp(m);
242#else
243 SymEigen< UpperSymmetrizeOperator<Matrix> > decomp(m.upperSymmetrize());
244#endif
245 if (!decomp.run()) return Type(0);
246 // compute (generalized) inverse matrix
247 decomp.ginv(inv);
248 return decomp.det();
249 }
250
251};
252
256template<class Matrix, int Size_>
257struct InvertUpperSymMatrixDispatcher
258{
259 typedef typename Matrix::Type Type;
271};
272
273template<class Matrix>
274struct InvertUpperSymMatrixDispatcher<Matrix, 1>
275{
276 typedef typename Matrix::Type Type;
279};
280
281template<class Matrix>
282struct InvertUpperSymMatrixDispatcher<Matrix, 2>
283{
284 typedef typename Matrix::Type Type;
287};
288
289template<class Matrix>
290struct InvertUpperSymMatrixDispatcher<Matrix, 3>
291{
292 typedef typename Matrix::Type Type;
295};
296
297template<class Matrix>
298struct InvertUpperSymMatrixDispatcher<Matrix, 4>
299{
300 typedef typename Matrix::Type Type;
303};
304
305} // namespace hidden
306
307} // namespace STK
308
309#endif /* STK_UPPERSYMINVERT_H */
STK_SymEigen.h
In this file we define the SymEigen class (for a symmetric matrix).
STK_lapack_SymEigen.h
In this file we define the enclosing class of the syevr lapck routine.
STK::ISymEigen::ginv
ArraySquare & ginv(ArraySquare &res) const
Compute the generalized inverse of the symmetric matrix and put the result in res.
Definition STK_ISymEigen.h:156
STK::ISymEigen::run
virtual bool run()
Find the eigenvalues and eigenvectors of the matrix.
Definition STK_ISymEigen.h:189
STK::ISymEigen::det
Type const & det() const
Definition STK_ISymEigen.h:127
STK::MultidimRegression
The MultidimRegression class allows to regress a multidimensional output variable among a multivariat...
Definition STK_MultidimRegression.h:52
STK::hidden::InvertUpperSymMatrixImpl::invertUpperSymMatrix22
static Type invertUpperSymMatrix22(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 2x2 matrix m using only its upper part and store the result in i...
Definition STK_InvertUpperSym.h:82
STK::hidden::InvertUpperSymMatrixImpl::invertUpperSymMatrix44
static Type invertUpperSymMatrix44(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 4x4 matrix m using only its upper part and store the result in i...
Definition STK_InvertUpperSym.h:143
STK::hidden::InvertUpperSymMatrixImpl::invertUpperSymMatrix11
static Type invertUpperSymMatrix11(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric matrix m of size 1x1 and store the result in inv.
Definition STK_InvertUpperSym.h:62
STK::hidden::InvertUpperSymMatrixImpl::invertUpperSymMatrixXX
static Type invertUpperSymMatrixXX(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 4x4 matrix m using only its upper part and store the result in i...
Definition STK_InvertUpperSym.h:238
STK::hidden::InvertUpperSymMatrixImpl::invertUpperSymMatrix33
static Type invertUpperSymMatrix33(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 3x3 matrix m using only its upper part and store the result in i...
Definition STK_InvertUpperSym.h:108
STK
The namespace STK is the main domain space of the Statistical ToolKit project.
STK::hidden::InvertUpperSymMatrixDispatcher< Matrix, 1 >::run
static Type run(Matrix const &m, CArraySquare< Type, 1 > &inv)
Definition STK_InvertUpperSym.h:277
STK::hidden::InvertUpperSymMatrixDispatcher< Matrix, 2 >::run
static Type run(Matrix const &m, CArraySquare< Type, 2 > &inv)
Definition STK_InvertUpperSym.h:285
STK::hidden::InvertUpperSymMatrixDispatcher< Matrix, 3 >::run
static Type run(Matrix const &m, CArraySquare< Type, 3 > &inv)
Definition STK_InvertUpperSym.h:293
STK::hidden::InvertUpperSymMatrixDispatcher< Matrix, 4 >::run
static Type run(Matrix const &m, CArraySquare< Type, 4 > &inv)
Definition STK_InvertUpperSym.h:301
STK::hidden::InvertUpperSymMatrixDispatcher
computing the inverse of a symmetric matrix.
Definition STK_InvertUpperSym.h:258
STK::hidden::InvertUpperSymMatrixDispatcher::run
static Type run(Matrix const &m, CArraySquare< Type, Size_ > &inv)
Definition STK_InvertUpperSym.h:260
STK::hidden::InvertUpperSymMatrixImpl
Implementation of the inversion matrix method for upper symmetric matrices.
Definition STK_InvertUpperSym.h:54