1 // -*- C++ -*-
2 /**
3 * @brief MatrixImplementation implements the classical mathematical Matrix
4 *
5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6 *
7 * This library is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library. If not, see <http://www.gnu.org/licenses/>.
19 *
20 */
21 #ifndef OPENTURNS_MATRIXIMPLEMENTATION_HXX
22 #define OPENTURNS_MATRIXIMPLEMENTATION_HXX
23
24 #include "openturns/PersistentCollection.hxx"
25 #include "openturns/Description.hxx"
26 #include "openturns/Point.hxx"
27 #include "openturns/Collection.hxx"
28
29 BEGIN_NAMESPACE_OPENTURNS
30
31 /**
32 * @class MatrixImplementation
33 *
34 * MatrixImplementation implements the classical mathematical MatrixImplementation
35 */
36
37 // Forward declaration of ComplexMatrixImplementation
38 class ComplexMatrixImplementation;
39 class Sample;
40
41 class OT_API MatrixImplementation
42 : public PersistentCollection<Scalar>
43
44 {
45 CLASSNAME
46
47 #ifndef SWIG
48 /** Declaration of friend operators */
operator *(const Scalar s,const MatrixImplementation & matrix)49 friend MatrixImplementation operator * (const Scalar s,
50 const MatrixImplementation & matrix)
51 {
52 return matrix.operator * (s);
53 }
54 #endif
55
56
57 public:
58
59 typedef Collection<Scalar> ScalarCollection;
60 typedef Collection<Complex> ComplexCollection;
61
62 /** Default constructor */
63 MatrixImplementation();
64
65 /** Constructor with size (rowDim and colDim) */
66 MatrixImplementation(const UnsignedInteger rowDim,
67 const UnsignedInteger colDim);
68
69 /** Constructor from range of external collection */
70 template <class InputIterator>
71 MatrixImplementation(const UnsignedInteger rowDim,
72 const UnsignedInteger colDim,
73 const InputIterator first,
74 const InputIterator last);
75
76 /** Constructor from external collection */
77 /** If the dimensions of the matrix and of the collection */
78 /** do not correspond, either the collection is truncated */
79 /** or the rest of the matrix is filled with zeros */
80 MatrixImplementation(const UnsignedInteger rowDim,
81 const UnsignedInteger colDim,
82 const ScalarCollection & elementsValues);
83
84 /** Virtual constructor */
85 MatrixImplementation * clone() const override;
86
87 /** String converter */
88 String __repr__() const override;
89 String __str__(const String & offset = "") const override;
90
91 /** Operator () gives access to the elements of the MatrixImplementation (to modify these elements) */
92 /** The element of the MatrixImplementation is designated by its row number i and its column number j */
93 Scalar & operator () (const UnsignedInteger i,
94 const UnsignedInteger j);
95
96 /** Operator () gives access to the elements of the MatrixImplementation (read only) */
97 /** The element of the MatrixImplementation is designated by its row number i and its column number j */
98 const Scalar & operator () (const UnsignedInteger i,
99 const UnsignedInteger j) const;
100
101 /** Get the dimensions of the MatrixImplementation */
102 /** Number of rows */
103 UnsignedInteger getNbRows() const;
104 /** Number of columns */
105 UnsignedInteger getNbColumns() const;
106 /** Dimension (for square matrices only */
107 UnsignedInteger getDimension() const;
108
109 /** MatrixImplementation transpose */
110 MatrixImplementation transpose() const;
111
112 /** MatrixImplementation reshape */
113 MatrixImplementation reshape(const UnsignedInteger newRowDim,
114 const UnsignedInteger newColDim) const;
115 void reshapeInPlace(const UnsignedInteger newRowDim,
116 const UnsignedInteger newColDim);
117
118
119 /** Row extraction */
120 const MatrixImplementation getRow(const UnsignedInteger rowIndex) const;
121 const MatrixImplementation getRowSym(const UnsignedInteger rowIndex) const;
122 /** Column extration */
123 const MatrixImplementation getColumn(const UnsignedInteger columnIndex) const;
124 const MatrixImplementation getColumnSym(const UnsignedInteger columnIndex) const;
125
126 /** MatrixImplementation addition (must have the same dimensions) */
127 MatrixImplementation operator + (const MatrixImplementation & matrix) const;
128
129 /** In-place MatrixImplementation addition (must have the same dimensions) */
130 MatrixImplementation & operator += (const MatrixImplementation & matrix);
131
132 /** MatrixImplementation subtraction (must have the same dimensions) */
133 MatrixImplementation operator - (const MatrixImplementation & matrix) const;
134
135 /** In-place MatrixImplementation subtraction (must have the same dimensions) */
136 MatrixImplementation & operator -= (const MatrixImplementation & matrix);
137
138 /** MatrixImplementation multiplications (must have consistent dimensions) */
139 MatrixImplementation genProd (const MatrixImplementation & matrix,
140 const Bool transposeLeft = false,
141 const Bool transposeRight = false) const;
142 MatrixImplementation symProd (const MatrixImplementation & m,
143 const char symSide) const;
144
145 /** MatrixImplementation integer power */
146 MatrixImplementation genPower(const UnsignedInteger n) const;
147 MatrixImplementation symPower(const UnsignedInteger n) const;
148
149 /** Multiplications with a SampleImplementation (must have consistent dimensions) */
150 Sample genSampleProd (const Sample & sample,
151 const Bool transposeMatrix,
152 const Bool transposeSample,
153 const char side) const;
154
155 /** Multiplications with a Point (must have consistent dimensions) */
156 Point genVectProd (const Point & pt,
157 const Bool transpose = false) const;
158 Point symVectProd (const Point & pt) const;
159
160 /** Using triangular matrix */
161 ScalarCollection triangularVectProd (const ScalarCollection & pt,
162 const char side = 'L',
163 const Bool transpose = false) const;
164 ScalarCollection triangularVectProd (const Point & pt,
165 const char side = 'L',
166 const Bool transpose = false) const;
167
168 /** Multiplication with a Scalar */
169 MatrixImplementation operator * (const Scalar s) const;
170
171 /** In-place Multiplication with a Scalar */
172 MatrixImplementation & operator *= (const Scalar s);
173
174 /** Division by a Scalar*/
175 MatrixImplementation operator / (const Scalar s) const;
176
177 /** In-place Division by a Scalar*/
178 MatrixImplementation & operator /= (const Scalar s);
179
180 /** Symmetrize MatrixImplementation in case it is a symmetric matrix (stored as a triangular matrix) */
181 void symmetrize() const;
182
183 /** Triangularize MatrixImplementation in case it is a triangular matrix (stored as a square matrix) */
184 void triangularize(const Bool isLowerTriangular) const;
185
186 /** Resolution of a linear system in case of a rectangular matrix */
187 Point solveLinearSystemRect(const Point & b,
188 const Bool keepIntact = true);
189 MatrixImplementation solveLinearSystemRect(const MatrixImplementation & b,
190 const Bool keepIntact = true);
191
192 /** Resolution of a linear system in case of a square matrix */
193 Point solveLinearSystemSquare(const Point & b,
194 const Bool keepIntact = true);
195 MatrixImplementation solveLinearSystemSquare(const MatrixImplementation & b,
196 const Bool keepIntact = true);
197
198 /** Resolution of a linear system in case of a triangular matrix */
199 Point solveLinearSystemTri(const Point & b,
200 const Bool keepIntact = true,
201 const Bool lower = true,
202 const Bool transpose = false);
203
204 MatrixImplementation solveLinearSystemTri(const MatrixImplementation & b,
205 const Bool keepIntact = true,
206 const Bool lower = true,
207 const Bool transpose = false);
208
209 /** Resolution of a linear system in case of a symmetric matrix */
210 Point solveLinearSystemSym(const Point & b,
211 const Bool keepIntact = true);
212 MatrixImplementation solveLinearSystemSym(const MatrixImplementation & b,
213 const Bool keepIntact = true);
214
215 /** Resolution of a linear system in case of a covariance matrix */
216 Point solveLinearSystemCov(const Point & b,
217 const Bool keepIntact = true);
218 MatrixImplementation solveLinearSystemCov(const MatrixImplementation & b,
219 const Bool keepIntact = true);
220
221 /** Triangular matrix product : side argument L/R for the position of the triangular matrix, up/lo to tell if it */
222 MatrixImplementation triangularProd(const MatrixImplementation & m,
223 const char side = 'L',
224 const char uplo = 'L') const;
225
226 /** Compute determinant */
227 Scalar computeLogAbsoluteDeterminant(Scalar & signOut,
228 const Bool keepIntact = true);
229 Scalar computeDeterminant(const Bool keepIntact = true);
230 Scalar computeLogAbsoluteDeterminantSym(Scalar & signOut,
231 const Bool keepIntact = true);
232 Scalar computeDeterminantSym(const Bool keepIntact = true);
233
234 /** Compute trace */
235 Scalar computeTrace() const;
236
237 /** Compute eigenvalues */
238 ComplexCollection computeEigenValuesSquare(const Bool keepIntact = true);
239 ComplexCollection computeEVSquare(ComplexMatrixImplementation & vOut,
240 const Bool keepIntact = true);
241 Point computeEigenValuesSym(const Bool keepIntact = true);
242 Point computeEVSym(MatrixImplementation & vOut,
243 const Bool keepIntact = true);
244 /** Compute the largest eigenvalue module using power iterations */
245 Bool computeLargestEigenValueModuleSquare(Scalar & maximumModule,
246 const UnsignedInteger maximumIterations = ResourceMap::GetAsUnsignedInteger("Matrix-LargestEigenValueIterations"),
247 const Scalar epsilon = ResourceMap::GetAsScalar("Matrix-LargestEigenValueRelativeError")) const;
248 Bool computeLargestEigenValueModuleSym(Scalar & maximumModule,
249 const UnsignedInteger maximumIterations = ResourceMap::GetAsUnsignedInteger("Matrix-LargestEigenValueIterations"),
250 const Scalar epsilon = ResourceMap::GetAsScalar("Matrix-LargestEigenValueRelativeError")) const;
251
252 /** Compute singular values */
253 Point computeSingularValues(const Bool keepIntact = true);
254
255 /** Build the singular value decomposition */
256 Point computeSVD(MatrixImplementation & uOut,
257 MatrixImplementation & vTOut,
258 const Bool fullSVD = false,
259 const Bool keepIntact = true);
260
261 /** Check if the matrix is symmetric */
262 virtual Bool isSymmetric() const;
263
264 /** Check if the matrix is SPD */
265 virtual Bool isPositiveDefinite() const;
266
267 /** Check if the matrix values belong to (-1;1) */
268 virtual Bool hasUnitRange() const;
269
270 /** Set small elements to zero */
271 virtual MatrixImplementation clean(const Scalar threshold) const;
272
273 virtual MatrixImplementation cleanSym(const Scalar threshold) const;
274
275 /** Build the Cholesky factorization of the matrix */
276 virtual MatrixImplementation computeCholesky(const Bool keepIntact = true);
277
278 #ifndef SWIG
279 /** Update in-place the Cholesky factor L of an SPD matrix M given a rank-one update vv^T, ie L becomes Lnew such that LnewLnew^t = Mnew with Mnew = M + vv^t */
280 static void CholeskyUpdate(MatrixImplementation & cholesky,
281 const Point & vector);
282
283 /** Downdate in-place the Cholesky factor L of an SPD matrix M given a rank-one downdate vv^T, ie L becomes Lnew such that LnewLnew^t = Mnew with Mnew = M - vv^t */
284 static void CholeskyDowndate(MatrixImplementation & cholesky,
285 const Point & vector);
286 #endif
287
288 /** Build the QR factorization of the matrix */
289 virtual MatrixImplementation computeQR(MatrixImplementation & ROut,
290 const Bool fullQR = false,
291 const Bool keepIntact = true);
292
293 /** Compute the Gram matrix associated to the matrix. If transpose == true, compute M^T.M, else M.M^T. */
294 virtual MatrixImplementation computeGram(const Bool transpose = true) const;
295
296 /** Comparison operators */
297 Bool operator == (const MatrixImplementation & rhs) const;
operator !=(const MatrixImplementation & rhs) const298 inline Bool operator != (const MatrixImplementation & rhs) const
299 {
300 return !((*this) == rhs);
301 }
302
303 /** Empty returns true if there is no element in the MatrixImplementation */
304 Bool isEmpty() const;
305
306 /** Returns true if triangular lower or upper */
307 Bool isTriangular(Bool lower = true) const;
308
309 /** Method save() stores the object through the StorageManager */
310 void save(Advocate & adv) const override;
311
312 /** Method load() reloads the object from the StorageManager */
313 void load(Advocate & adv) override;
314
315 /** Low-level data access */
316 UnsignedInteger stride(const UnsignedInteger dim) const;
317
318 protected:
319
320 /** MatrixImplementation Dimensions */
321 UnsignedInteger nbRows_;
322 UnsignedInteger nbColumns_;
323
324 /** Position conversion function : the indices i & j are used to compute the actual position of the element in the collection */
325 inline UnsignedInteger convertPosition (const UnsignedInteger i,
326 const UnsignedInteger j) const;
327
328
329
330 }; /* class MatrixImplementation */
331
convertPosition(const UnsignedInteger i,const UnsignedInteger j) const332 inline UnsignedInteger MatrixImplementation::convertPosition (const UnsignedInteger i,
333 const UnsignedInteger j) const
334 {
335 return i + nbRows_ * j ;
336 }
337
338 /** Constructor from range of external collection */
339 template <class InputIterator>
MatrixImplementation(const UnsignedInteger rowDim,const UnsignedInteger colDim,const InputIterator first,const InputIterator last)340 MatrixImplementation::MatrixImplementation(const UnsignedInteger rowDim,
341 const UnsignedInteger colDim,
342 const InputIterator first,
343 const InputIterator last)
344 : PersistentCollection<Scalar>(rowDim * colDim, 0.0),
345 nbRows_(rowDim),
346 nbColumns_(colDim)
347 {
348 this->assign(first, last);
349 }
350
351 END_NAMESPACE_OPENTURNS
352
353 #endif /* OPENTURNS_MATRIXIMPLEMENTATION_HXX */
354