STK++ 0.9.13
STK_InvertLowerSym.h
Go to the documentation of this file.
1/*--------------------------------------------------------------------*/
2/* Copyright (C) 2004-2016 Serge Iovleff, jShift+Université Lille 1, jShift+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, jShift+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, jShift+write to the
16 Free Software Foundation, jShift+Inc.,
17 59 Temple Place,
18 Suite 330,
19 Boston, jShift+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, jShift+S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
29 **/
30
35#ifndef STK_INVERTLOWERSYM_H
36#define STK_INVERTLOWERSYM_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 InvertLowerSymMatrixImpl
54{
55 typedef typename Matrix::Type Type;
62 static Type invertLowerSymMatrix11( 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 invertLowerSymMatrix22( Matrix const& m, CArraySquare<Type, Size_>& inv)
83 {
84 const int iShift = m.beginRows(), jShift = m.beginCols();
85 inv.resize(TRange<Size_>(0, 1));
86 // cofactor (0,0) [0]
87 inv(0, 0) = m(iShift+1, jShift+1);
88 // cofactor (1,0) [1]
89 inv(1, 0) = - m(iShift+1, jShift+0);
90 // cofactor (1,1) [3]
91 inv(1, 1) = m(iShift+0, jShift+0);
92 // symmetry
93 inv(0, 1) = inv(1, 0);
94 // determinant
95 Type det = m(iShift+0, jShift+0) * inv(0, 0)
96 + m(iShift+1, jShift+0) * inv(1, 0);
97 if (det == Type(0)) return Type(0);
98 // inverse matrix
99 inv /= det;
100 return det;
101 }
108 static Type invertLowerSymMatrix33( Matrix const& m, CArraySquare<Type, Size_>& inv)
109 {
110 const int iShift = m.beginRows(), jShift = m.beginCols();
111 inv.resize(TRange<Size_>(0, 2));
112 // cofactor (0,0) [0]
113 inv(0, 0) = m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
114 - m(iShift+2, jShift+1) * m(iShift+2, jShift+1);
115 // cofactor (1,0) [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(2, 0) = m(iShift+1, jShift+0) * m(iShift+2, jShift+1)
119 - m(iShift+2, jShift+0) * 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+2, jShift+0);
122 inv(2, 1) = m(iShift+2, jShift+0) * m(iShift+1, jShift+0)
123 - m(iShift+0, jShift+0) * m(iShift+2, jShift+1);
124 inv(2, 2) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1)
125 - m(iShift+1, jShift+0) * m(iShift+0, jShift+1);
126 // symmetry
127 inv(0,1) = inv(1,0); inv(0,2) = inv(2,0);
128 inv(1,2) = inv(2,1);
129
130 // compute determinant and inverse matrix
131 Type det = m(iShift+0, jShift+0) * inv(0, 0)
132 + m(iShift+1, jShift+0) * inv(1, 0)
133 + m(iShift+2, jShift+0) * inv(2,0);
134 if (det == Type(0)) return Type(0);
135 // inverse matrix
136 inv /= det;
137 return det;
138 }
145 static Type invertLowerSymMatrix44( Matrix const& m, CArraySquare<Type, Size_>& inv)
146 {
147 const int iShift = m.beginRows(), jShift = m.beginCols();
148 inv.resize(TRange<Size_>(0, jShift+3));
149 // cofactor (0,0) [0]
150 inv(0,0) = m(iShift+1, jShift+1) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
151 - m(iShift+1, jShift+1) * m(iShift+3, jShift+2) * m(iShift+3, jShift+2)
152 - m(iShift+2, jShift+1) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
153 + m(iShift+2, jShift+1) * m(iShift+3, jShift+1) * m(iShift+3, jShift+2)
154 + m(iShift+3, jShift+1) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
155 - m(iShift+3, jShift+1) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2);
156 // cofactor (1,0) [1]
157 inv(1,0) = -m(iShift+1, jShift+0) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
158 + m(iShift+1, jShift+0) * m(iShift+3, jShift+2) * m(iShift+3, jShift+2)
159 + m(iShift+2, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
160 - m(iShift+2, jShift+1) * m(iShift+3, jShift+0) * m(iShift+3, jShift+2)
161 - m(iShift+3, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
162 + m(iShift+3, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2);
163 // cofactor (2,0) [2]
164 inv(2,0) = m(iShift+1, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
165 - m(iShift+1, jShift+0) * m(iShift+3, jShift+1) * m(iShift+3, jShift+2)
166 - m(iShift+1, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
167 + m(iShift+1, jShift+1) * m(iShift+3, jShift+0) * m(iShift+3, jShift+2)
168 + m(iShift+3, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
169 - m(iShift+3, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
170 // cofactor (3,0) [3]
171 inv(3,0) = -m(iShift+1, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
172 + m(iShift+1, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2)
173 + m(iShift+1, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
174 - m(iShift+1, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2)
175 - m(iShift+2, jShift+1) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
176 + m(iShift+2, jShift+1) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
177 // cofactor (1,1) [5]
178 inv(1,1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+2) * m(iShift+3, jShift+3)
179 - m(iShift+0, jShift+0) * m(iShift+3, jShift+2) * m(iShift+3, jShift+2)
180 - m(iShift+2, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
181 + m(iShift+2, jShift+0) * m(iShift+3, jShift+0) * m(iShift+3, jShift+2)
182 + m(iShift+3, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
183 - m(iShift+3, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2);
184 // cofactor (2,1) [6]
185 inv(2,1) = -m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+3)
186 + m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+3, jShift+2)
187 + m(iShift+1, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+3)
188 - m(iShift+1, jShift+0) * m(iShift+3, jShift+0) * m(iShift+3, jShift+2)
189 - m(iShift+3, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
190 + m(iShift+3, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
191 // cofactor (3,1) [7]
192 inv(3,1) = m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+3, jShift+2)
193 - m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+2)
194 - m(iShift+1, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+2)
195 + m(iShift+1, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+2)
196 + m(iShift+2, jShift+0) * m(iShift+2, jShift+0) * m(iShift+3, jShift+1)
197 - m(iShift+2, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1);
198 // cofactor (2,2) [10]
199 inv(2,2) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+3, jShift+3)
200 - m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+3, jShift+1)
201 - m(iShift+1, jShift+0) * m(iShift+1, jShift+0) * m(iShift+3, jShift+3)
202 + m(iShift+1, jShift+0) * m(iShift+3, jShift+0) * m(iShift+3, jShift+1)
203 + m(iShift+3, jShift+0) * m(iShift+1, jShift+0) * m(iShift+3, jShift+1)
204 - m(iShift+3, jShift+0) * m(iShift+3, jShift+0) * m(iShift+1, jShift+1);
205 // cofactor (3,2) [11]
206 inv(3,2) = -m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+3, jShift+2)
207 + m(iShift+0, jShift+0) * m(iShift+3, jShift+1) * m(iShift+2, jShift+1)
208 + m(iShift+1, jShift+0) * m(iShift+1, jShift+0) * m(iShift+3, jShift+2)
209 - m(iShift+1, jShift+0) * m(iShift+3, jShift+0) * m(iShift+2, jShift+1)
210 - m(iShift+2, jShift+0) * m(iShift+1, jShift+0) * m(iShift+3, jShift+1)
211 + m(iShift+2, jShift+0) * m(iShift+3, jShift+0) * m(iShift+1, jShift+1);
212 // cofactor (3,3) [15]
213 inv(3,3) = m(iShift+0, jShift+0) * m(iShift+1, jShift+1) * m(iShift+2, jShift+2)
214 - m(iShift+0, jShift+0) * m(iShift+2, jShift+1) * m(iShift+2, jShift+1)
215 - m(iShift+1, jShift+0) * m(iShift+1, jShift+0) * m(iShift+2, jShift+2)
216 + m(iShift+1, jShift+0) * m(iShift+2, jShift+0) * m(iShift+2, jShift+1)
217 + m(iShift+2, jShift+0) * m(iShift+1, jShift+0) * m(iShift+2, jShift+1)
218 - m(iShift+2, jShift+0) * m(iShift+2, jShift+0) * m(iShift+1, jShift+1);
219 // symmetric
220 inv(0,1) = inv(1,0); inv(0,2) = inv(2,0); inv(0,3) = inv(3,0);
221 inv(1,2) = inv(2,1); inv(1,3) = inv(3,1);
222 inv(2,3) = inv(3,2);
223 // compute determinant
224 Type det = m(iShift+0, jShift+0) * inv(0,0)
225 + m(iShift+1, jShift+0) * inv(1,0)
226 + m(iShift+2, jShift+0) * inv(2,0)
227 + m(iShift+3, jShift+0) * inv(3,0);
228 if (det == Type(0)) return Type(0);
229 // inverse matrix
230 inv /= det;
231 return det;
232 }
233
241 static Type invertLowerSymMatrixXX( Matrix const& m, CArraySquare<Type, Size_>& inv)
242 {
243#ifdef STKUSELAPACK
244 lapack::SymEigen<Matrix> decomp(m);
245 decomp.setUplo('L'); // default is U
246#else
247 SymEigen< LowerSymmetrizeOperator<Matrix> > decomp(m.lowerSymmetrize());
248#endif
249 if (!decomp.run()) return Type(0);
250 decomp.ginv(inv);
251 return decomp.det();
252 }
253};
254
258template<class Matrix, int Size_>
259struct InvertLowerSymMatrixDispatcher
260{
261 typedef typename Matrix::Type Type;
273};
274
275template<class Matrix>
276struct InvertLowerSymMatrixDispatcher<Matrix, 1>
277{
278 typedef typename Matrix::Type Type;
281};
282
283template<class Matrix>
284struct InvertLowerSymMatrixDispatcher<Matrix, 2>
285{
286 typedef typename Matrix::Type Type;
289};
290
291template<class Matrix>
292struct InvertLowerSymMatrixDispatcher<Matrix, 3>
293{
294 typedef typename Matrix::Type Type;
297};
298
299template<class Matrix>
300struct InvertLowerSymMatrixDispatcher<Matrix, 4>
301{
302 typedef typename Matrix::Type Type;
305};
306
307} // namespace hidden
308
309} // namespace STK
310
311#endif /* STK_INVERTLOWERSYM_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::InvertLowerSymMatrixImpl::invertLowerSymMatrix33
static Type invertLowerSymMatrix33(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 3x3 matrix m using only its lower part and store the result in i...
Definition STK_InvertLowerSym.h:108
STK::hidden::InvertLowerSymMatrixImpl::invertLowerSymMatrix22
static Type invertLowerSymMatrix22(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the lower symmetric 2x2 matrix m using only its lower part and store the resul...
Definition STK_InvertLowerSym.h:82
STK::hidden::InvertLowerSymMatrixImpl::invertLowerSymMatrix44
static Type invertLowerSymMatrix44(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the symmetric 4x4 matrix m using only its lower part and store the result in i...
Definition STK_InvertLowerSym.h:145
STK::hidden::InvertLowerSymMatrixImpl::invertLowerSymMatrixXX
static Type invertLowerSymMatrixXX(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the lower symmetric 4x4 matrix m using only its lower part and store the resul...
Definition STK_InvertLowerSym.h:241
STK::hidden::InvertLowerSymMatrixImpl::invertLowerSymMatrix11
static Type invertLowerSymMatrix11(Matrix const &m, CArraySquare< Type, Size_ > &inv)
compute the inverse of the lower symmetric matrix m of size 1x1 and store the result in inv.
Definition STK_InvertLowerSym.h:62
STK
The namespace STK is the main domain space of the Statistical ToolKit project.
STK::hidden::InvertLowerSymMatrixDispatcher< Matrix, 1 >::run
static Type run(Matrix const &m, CArraySquare< Type, 1 > &inv)
Definition STK_InvertLowerSym.h:279
STK::hidden::InvertLowerSymMatrixDispatcher< Matrix, 2 >::run
static Type run(Matrix const &m, CArraySquare< Type, 2 > &inv)
Definition STK_InvertLowerSym.h:287
STK::hidden::InvertLowerSymMatrixDispatcher< Matrix, 3 >::run
static Type run(Matrix const &m, CArraySquare< Type, 3 > &inv)
Definition STK_InvertLowerSym.h:295
STK::hidden::InvertLowerSymMatrixDispatcher< Matrix, 4 >::run
static Type run(Matrix const &m, CArraySquare< Type, 4 > &inv)
Definition STK_InvertLowerSym.h:303
STK::hidden::InvertLowerSymMatrixDispatcher
computing the inverse of a lower symmetric matrix.
Definition STK_InvertLowerSym.h:260
STK::hidden::InvertLowerSymMatrixDispatcher::run
static Type run(Matrix const &m, CArraySquare< Type, Size_ > &inv)
Definition STK_InvertLowerSym.h:262
STK::hidden::InvertLowerSymMatrixImpl
Implementation of the computation of the inverse of a lower symmetric matrix.
Definition STK_InvertLowerSym.h:54