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