1 // -*- C++ -*-
2 /**
3 * @brief MatrixImplementation implements the classical mathematical MatrixImplementation
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 #include <cstdlib>
22 #include <functional>
23
24 #include "openturns/MatrixImplementation.hxx"
25 #include "openturns/ComplexMatrixImplementation.hxx"
26 #include "openturns/PersistentObjectFactory.hxx"
27 #include "openturns/Lapack.hxx"
28 #include "openturns/ResourceMap.hxx"
29 #include "openturns/SpecFunc.hxx"
30 #include "openturns/Os.hxx"
31 #include "openturns/Exception.hxx"
32 #include "openturns/Sample.hxx"
33
34 BEGIN_NAMESPACE_OPENTURNS
35
36 CLASSNAMEINIT(MatrixImplementation)
37
38
39 static const Factory<MatrixImplementation> Factory_MatrixImplementation;
40
41 // All the pivots with a magnitude less than this threshold are considered as zero
42 /* Default constructor */
MatrixImplementation()43 MatrixImplementation::MatrixImplementation()
44 : PersistentCollection<Scalar>()
45 , nbRows_(0)
46 , nbColumns_(0)
47 {
48 // Nothing to do
49 }
50
51 /* Constructor with size (rowDim and colDim) */
52 /* The MatrixImplementation is made up of a collection of rowDim*colDim elements */
53 /* The MatrixImplementation is viewed as a set of column vectors read one after another */
MatrixImplementation(const UnsignedInteger rowDim,const UnsignedInteger colDim)54 MatrixImplementation::MatrixImplementation(const UnsignedInteger rowDim,
55 const UnsignedInteger colDim)
56 : PersistentCollection<Scalar>(rowDim * colDim, 0.0)
57 , nbRows_(rowDim)
58 , nbColumns_(colDim)
59 {
60 // Nothing to do
61 }
62
63 /* Constructor from external collection */
MatrixImplementation(const UnsignedInteger rowDim,const UnsignedInteger colDim,const Collection<Scalar> & elementsValues)64 MatrixImplementation::MatrixImplementation(const UnsignedInteger rowDim,
65 const UnsignedInteger colDim,
66 const Collection<Scalar> & elementsValues)
67 : PersistentCollection<Scalar>(rowDim * colDim, 0.0)
68 , nbRows_(rowDim)
69 , nbColumns_(colDim)
70 {
71 const UnsignedInteger matrixSize = std::min(rowDim * colDim, elementsValues.getSize());
72 std::copy(elementsValues.begin(), elementsValues.begin() + matrixSize, begin());
73 }
74
75 /* Virtual constructor */
clone() const76 MatrixImplementation * MatrixImplementation::clone() const
77 {
78 return new MatrixImplementation(*this);
79 }
80
81
82 /* String converter */
__repr__() const83 String MatrixImplementation::__repr__() const
84 {
85 return OSS(true) << "class=" << getClassName()
86 << " name=" << getName()
87 << " rows=" << nbRows_
88 << " columns=" << nbColumns_
89 << " values=" << PersistentCollection<Scalar>::__repr__();
90 }
91
__str__(const String & offset) const92 String MatrixImplementation::__str__(const String & offset) const
93 {
94 OSS oss(false);
95 if (nbRows_ == 0 || nbColumns_ == 0) return "[]";
96
97 if ( (nbRows_ >= ResourceMap::GetAsUnsignedInteger("Matrix-size-visible-in-str-from")) ||
98 (nbColumns_ >= ResourceMap::GetAsUnsignedInteger("Matrix-size-visible-in-str-from")) )
99 oss << nbRows_ << "x" << nbColumns_ << Os::GetEndOfLine();
100
101 size_t lwidth(0);
102 size_t rwidth(0);
103 for ( UnsignedInteger i = 0; i < nbRows_; ++i )
104 for ( UnsignedInteger j = 0; j < nbColumns_; ++j )
105 {
106 String st(OSS() << (*this)(i, j));
107 size_t dotpos(st.find( '.' ));
108 lwidth = std::max( lwidth, (dotpos != String::npos) ? dotpos : st.size() );
109 rwidth = std::max( rwidth, (dotpos != String::npos) ? st.size() - dotpos : 0 );
110 }
111
112 const char * bracket("[");
113 const char * newline("");
114 const char * noffset("");
115 for ( UnsignedInteger i = 0; i < nbRows_; ++i, newline = Os::GetEndOfLine(), noffset = offset.c_str(), bracket = " " )
116 {
117 oss << newline << noffset << bracket << "[ ";
118 const char * sep("");
119 for ( UnsignedInteger j = 0; j < nbColumns_; ++j, sep = " " )
120 {
121 String st(OSS() << (*this)(i, j));
122 size_t dotpos(st.find( '.' ));
123 oss << sep << String( lwidth - ((dotpos != String::npos) ? dotpos : st.size()), ' ' )
124 << st
125 << String( rwidth - ((dotpos != String::npos) ? st.size() - dotpos : 0), ' ' );
126 }
127 oss << " ]";
128 }
129 oss << "]";
130 return oss;
131 }
132
133 /* Operator () gives access to the elements of the MatrixImplementation (to modify these elements) */
134 /* The element of the MatrixImplementation is designated by its row number i and its column number j */
135 /* the first element of the MatrixImplementation is m(0,0) */
operator ()(const UnsignedInteger i,const UnsignedInteger j)136 Scalar & MatrixImplementation::operator() (const UnsignedInteger i,
137 const UnsignedInteger j)
138 {
139 if (!(i < nbRows_)) throw OutOfBoundException(HERE) << "i (" << i << ") must be less than row dim (" << nbRows_ << ")";
140 if (!(j < nbColumns_)) throw OutOfBoundException(HERE) << "j (" << j << ") must be less than column dim (" << nbColumns_ << ")";
141 return operator[](convertPosition(i, j));
142 }
143
144 /* Operator () gives access to the elements of the MatrixImplementation (read only) */
145 /* The element of the MatrixImplementation is designated by its row number i and its column number j */
operator ()(const UnsignedInteger i,const UnsignedInteger j) const146 const Scalar & MatrixImplementation::operator() (const UnsignedInteger i,
147 const UnsignedInteger j) const
148 {
149 if (!(i < nbRows_)) throw OutOfBoundException(HERE) << "i (" << i << ") must be less than row dim (" << nbRows_ << ")";
150 if (!(j < nbColumns_)) throw OutOfBoundException(HERE) << "j (" << j << ") must be less than column dim (" << nbColumns_ << ")";
151 return operator[](convertPosition(i, j));
152 }
153
154
155 /* Get the dimensions of the MatrixImplementation : number of rows */
getNbRows() const156 UnsignedInteger MatrixImplementation::getNbRows() const
157 {
158 return nbRows_;
159 }
160
161 /* Get the dimensions of the MatrixImplementation : number of columns */
getNbColumns() const162 UnsignedInteger MatrixImplementation::getNbColumns() const
163 {
164 return nbColumns_;
165 }
166
167 /* Get the dimensions of the MatrixImplementation : dimension (square matrix : nbRows_) */
getDimension() const168 UnsignedInteger MatrixImplementation::getDimension() const
169 {
170 return nbRows_;
171 }
172
173 /* MatrixImplementation transpose */
transpose() const174 MatrixImplementation MatrixImplementation::transpose() const
175 {
176 MatrixImplementation trans(nbColumns_, nbRows_);
177 // The source matrix is accessed columnwise in the natural order
178 for (UnsignedInteger j = 0; j < nbColumns_; ++j)
179 for (UnsignedInteger i = 0; i < nbRows_; ++i)
180 trans(j, i) = operator()(i, j);
181 return trans;
182 }
183
184 /* MatrixImplementation reshape */
reshape(const UnsignedInteger newRowDim,const UnsignedInteger newColDim) const185 MatrixImplementation MatrixImplementation::reshape(const UnsignedInteger newRowDim,
186 const UnsignedInteger newColDim) const
187 {
188 return MatrixImplementation(newRowDim, newColDim, *this);
189 }
190
reshapeInPlace(const UnsignedInteger newRowDim,const UnsignedInteger newColDim)191 void MatrixImplementation::reshapeInPlace(const UnsignedInteger newRowDim,
192 const UnsignedInteger newColDim)
193 {
194 if (newRowDim * newColDim != getSize()) resize(newRowDim * newColDim);
195 nbRows_ = newRowDim;
196 nbColumns_ = newColDim;
197 }
198
199
200 /* Row extraction */
getRow(const UnsignedInteger rowIndex) const201 const MatrixImplementation MatrixImplementation::getRow(const UnsignedInteger rowIndex) const
202 {
203 if (!(rowIndex < nbRows_)) throw OutOfBoundException(HERE) << "Error: the row index=" << rowIndex << " must be less than the row number=" << nbRows_;
204 MatrixImplementation row(1, nbColumns_);
205 for (UnsignedInteger i = 0; i < nbColumns_; ++i) row(0, i) = (*this)(rowIndex, i);
206 return row;
207 }
208
getRowSym(const UnsignedInteger rowIndex) const209 const MatrixImplementation MatrixImplementation::getRowSym(const UnsignedInteger rowIndex) const
210 {
211 if (!(rowIndex < nbRows_)) throw OutOfBoundException(HERE) << "Error: the row index=" << rowIndex << " must be less than the row number=" << nbRows_;
212 MatrixImplementation row(1, nbColumns_);
213 for (UnsignedInteger i = 0; i < rowIndex; ++i) row(0, i) = (*this)(rowIndex, i);
214 for (UnsignedInteger i = rowIndex; i < nbColumns_; ++i) row(0, i) = (*this)(i, rowIndex);
215 return row;
216 }
217
218 /* Column extration */
getColumn(const UnsignedInteger columnIndex) const219 const MatrixImplementation MatrixImplementation::getColumn(const UnsignedInteger columnIndex) const
220 {
221 if (!(columnIndex < nbColumns_)) throw OutOfBoundException(HERE) << "Error: the column index=" << columnIndex << " must be less than the column number=" << nbColumns_;
222 MatrixImplementation column(nbRows_, 1);
223 for (UnsignedInteger i = 0; i < nbRows_; ++i) column(i, 0) = (*this)(i, columnIndex);
224 return column;
225 }
226
getColumnSym(const UnsignedInteger columnIndex) const227 const MatrixImplementation MatrixImplementation::getColumnSym(const UnsignedInteger columnIndex) const
228 {
229 if (!(columnIndex < nbColumns_)) throw OutOfBoundException(HERE) << "Error: the column index=" << columnIndex << " must be less than the column number=" << nbColumns_;
230 MatrixImplementation column(nbRows_, 1);
231 for (UnsignedInteger i = 0; i < columnIndex; ++i) column(i, 0) = (*this)(columnIndex, i);
232 for (UnsignedInteger i = columnIndex; i < nbRows_; ++i) column(i, 0) = (*this)(i, columnIndex);
233 return column;
234 }
235
236
237 /* MatrixImplementation addition (must have the same dimensions) */
operator +(const MatrixImplementation & matrix) const238 MatrixImplementation MatrixImplementation::operator+ (const MatrixImplementation & matrix) const
239 {
240 if ((nbRows_ != matrix.nbRows_ ) || (nbColumns_ != matrix.nbColumns_ )) throw InvalidDimensionException(HERE) << "Cannot add matrices with incompatible dimensions";
241 // Must copy as add will be overwritten by the operation
242 MatrixImplementation result(matrix);
243 int size(nbRows_ * nbColumns_);
244 double alpha(1.0);
245 int one(1);
246 daxpy_(&size, &alpha, const_cast<double*>(&((*this)[0])), &one, &result[0], &one);
247
248 return result;
249 }
250
251 /* In-place MatrixImplementation addition (must have the same dimensions) */
operator +=(const MatrixImplementation & matrix)252 MatrixImplementation & MatrixImplementation::operator+= (const MatrixImplementation & matrix)
253 {
254 if ((nbRows_ != matrix.nbRows_ ) || (nbColumns_ != matrix.nbColumns_ )) throw InvalidDimensionException(HERE) << "Cannot add matrices with incompatible dimensions";
255 // Must copy as add will be overwritten by the operation
256 MatrixImplementation result(matrix);
257 int size(nbRows_ * nbColumns_);
258 double alpha(1.0);
259 int one(1);
260 daxpy_(&size, &alpha, const_cast<double*>(&(matrix[0])), &one, &(*this)[0], &one);
261
262 return *this;
263 }
264
265 /* MatrixImplementation subtraction (must have the same dimensions) */
operator -(const MatrixImplementation & matrix) const266 MatrixImplementation MatrixImplementation::operator- (const MatrixImplementation & matrix) const
267 {
268 if ((nbRows_ != matrix.nbRows_ ) || (nbColumns_ != matrix.nbColumns_ )) throw InvalidDimensionException(HERE) << "Cannot subtract matrices with incompatible dimensions";
269 // Must copy as add will be overwritten by the operation
270 MatrixImplementation result(*this);
271 int size(nbRows_ * nbColumns_);
272 double alpha(-1.0);
273 int one(1);
274 daxpy_(&size, &alpha, const_cast<double*>(&(matrix[0])), &one, &result[0], &one);
275
276 return result;
277 }
278
279 /* In-place MatrixImplementation subtraction (must have the same dimensions) */
operator -=(const MatrixImplementation & matrix)280 MatrixImplementation & MatrixImplementation::operator-= (const MatrixImplementation & matrix)
281 {
282 if ((nbRows_ != matrix.nbRows_ ) || (nbColumns_ != matrix.nbColumns_ )) throw InvalidDimensionException(HERE) << "Cannot subtract matrices with incompatible dimensions";
283 // Must copy as add will be overwritten by the operation
284 MatrixImplementation result(*this);
285 int size(nbRows_ * nbColumns_);
286 double alpha(-1.0);
287 int one(1);
288 daxpy_(&size, &alpha, const_cast<double*>(&(matrix[0])), &one, &(*this)[0], &one);
289
290 return *this;
291 }
292
293 /* MatrixImplementation multiplications (must have consistent dimensions) */
genProd(const MatrixImplementation & matrix,const Bool transposeLeft,const Bool transposeRight) const294 MatrixImplementation MatrixImplementation::genProd (const MatrixImplementation & matrix,
295 const Bool transposeLeft,
296 const Bool transposeRight) const
297 {
298 // Multiply an mxk matrix op(A) by a kxn matrix op(B),
299 // where op(A) == A if transposeLeft == false, op(A) == A^t else
300 // and op(B) == B if transposeRight == false, op(B) == B^t else
301 int m(transposeLeft ? nbColumns_ : nbRows_);
302 int k(transposeLeft ? nbRows_ : nbColumns_);
303 int l(transposeRight ? matrix.nbColumns_ : matrix.nbRows_);
304 int n(transposeRight ? matrix.nbRows_ : matrix.nbColumns_);
305 if (k != l) throw InvalidDimensionException(HERE) << "Invalid dimensions in matrix/matrix product left="
306 << nbRows_ << "x" << nbColumns_
307 << " right=" << matrix.nbRows_ << "x"
308 << matrix.nbColumns_ << ", left is transposed=" << (transposeLeft ? "true" : "false") << ", right is transposed=" << (transposeRight ? "true" : "false");
309
310 MatrixImplementation mult(m, n);
311 if ((m == 0) || (n == 0) || (k == 0)) return mult;
312 char transa(transposeLeft ? 'T' : 'N');
313 char transb(transposeRight ? 'T' : 'N');
314 double alpha(1.0);
315 double beta(0.0);
316 int ltransa(1);
317 int ltransb(1);
318 int lda(nbRows_);
319 int ldb(matrix.nbRows_);
320 dgemm_(&transa, &transb, &m, &n, &k, &alpha, const_cast<double*>(&((*this)[0])), &lda, const_cast<double*>(&(matrix[0])), &ldb, &beta, &mult[0], &m, <ransa, <ransb);
321
322 return mult;
323 }
324
symProd(const MatrixImplementation & matrix,const char symSide) const325 MatrixImplementation MatrixImplementation::symProd (const MatrixImplementation & matrix,
326 const char symSide) const
327 {
328 const MatrixImplementation & left = (symSide == 'L') ? *this : matrix;
329 const MatrixImplementation & right = (symSide == 'L') ? matrix : *this;
330
331 if (left.nbColumns_ != right.nbRows_) throw InvalidDimensionException(HERE) << "Invalid dimensions in matrix/matrix product"
332 << " left=" << left.nbRows_ << "x" << left.nbColumns_
333 << " right=" << right.nbRows_ << "x" << right.nbColumns_;
334
335 MatrixImplementation mult(left.nbRows_, right.nbColumns_);
336 if ((left.nbRows_ == 0) || (left.nbColumns_ == 0) || (right.nbRows_ == 0) || (right.nbColumns_ == 0)) return mult;
337
338 char side = symSide;
339 char uplo = 'L';
340 int m = left.nbRows_;
341 int n = right.nbColumns_;
342 int lda = nbRows_;
343 int ldb = matrix.nbRows_;
344 double alpha = 1.0;
345 double beta = 0.0;
346 int lside = 1;
347 int luplo = 1;
348 dsymm_(&side, &uplo, &m, &n, &alpha, const_cast<double*>(&((*this)[0])), &lda, const_cast<double*>(&(matrix[0])), &ldb, &beta, &mult[0], &m, &lside, &luplo);
349 return mult;
350 }
351
352 /* Multiplications with a Point (must have consistent dimensions) */
genVectProd(const Point & pt,const Bool transposed) const353 Point MatrixImplementation::genVectProd (const Point & pt, const Bool transposed) const
354 {
355 int k = transposed ? nbRows_ : nbColumns_;
356 if (k != (int)pt.getDimension()) throw InvalidDimensionException(HERE) << "Invalid dimension in matrix/vector product: columns=" << k << " / vector dimension=" << pt.getDimension() << ".";
357
358 int l = transposed ? nbColumns_ : nbRows_;
359 Point prod(l);
360 if ((nbRows_ == 0) || (nbColumns_ == 0)) return prod;
361 char trans = transposed ? 'T' : 'N';
362 int one = 1;
363 double alpha = 1.0;
364 double beta = 0.0;
365 int ltrans = 1;
366 int m = nbRows_;
367 int n = nbColumns_;
368
369 dgemv_(&trans, &m, &n, &alpha, const_cast<double*>(&((*this)[0])), &m, const_cast<double*>(&(pt[0])), &one, &beta, &prod[0], &one, <rans);
370
371 return prod;
372 }
373
374 /* Multiplications with a SampleImplementation (must have consistent dimensions)
375 If side=='L', we compute this * S, otherwise we compute S * this */
genSampleProd(const Sample & sample,const Bool transposeMatrix,const Bool transposeSample,const char side) const376 Sample MatrixImplementation::genSampleProd (const Sample & sample, const Bool transposeMatrix, const Bool transposeSample, const char side) const
377 {
378 // Sample and matrices are stored in reverse order, thus instead of computing
379 // P = this * S, we compute P^T = S^T * this^T and S^T has a column-major order
380 // as expected by dgemm.
381 // Likewise when side is 'R', we compute P^T = this^T * S^T instead of P = S * this.
382 const UnsignedInteger matrixRows(transposeMatrix ? nbColumns_ : nbRows_);
383 const UnsignedInteger matrixColumns(transposeMatrix ? nbRows_ : nbColumns_);
384 const UnsignedInteger sampleRows(transposeSample ? sample.getDimension() : sample.getSize());
385 const UnsignedInteger sampleColumns(transposeSample ? sample.getSize() : sample.getDimension());
386 if (side == 'L' && matrixColumns != sampleRows) throw InvalidDimensionException(HERE) << "Invalid dimension in matrix*sample product: columns=" << matrixColumns << " / sample rows=" << sampleRows << ".";
387 if (side == 'R' && matrixRows != sampleColumns) throw InvalidDimensionException(HERE) << "Invalid dimension in sample*matrix product: rows=" << matrixRows << " / sample columns=" << sampleColumns << ".";
388
389 const UnsignedInteger prodRows(side == 'L' ? matrixRows : sampleRows);
390 const UnsignedInteger prodColumns(side == 'L' ? sampleColumns : matrixColumns);
391 Sample prod(prodRows, prodColumns);
392 if (nbRows_ == 0 || nbColumns_ == 0) return prod;
393 char transa = transposeSample ? 'T' : 'N';
394 char transb = transposeMatrix ? 'N' : 'T';
395 double alpha = 1.0;
396 double beta = 0.0;
397 int ltransa = 1;
398 int ltransb = 1;
399 if (side == 'L')
400 {
401 int m = sampleColumns;
402 int n = matrixRows;
403 int k = matrixColumns;
404 int lda = (transa == 'N' ? m : k);
405 int ldb = (transb == 'N' ? k : n);
406 dgemm_(&transa, &transb, &m, &n, &k, &alpha, const_cast<double*>(&(sample(0, 0))), &lda, const_cast<double*>(&((*this)[0])), &ldb, &beta, &prod(0, 0), &m, <ransa, <ransb);
407 }
408 else
409 {
410 int m = matrixColumns;
411 int n = sampleRows;
412 int k = matrixRows;
413 int lda = (transb == 'N' ? m : k);
414 int ldb = (transa == 'N' ? k : n);
415 dgemm_(&transb, &transa, &m, &n, &k, &alpha, const_cast<double*>(&((*this)[0])), &lda, const_cast<double*>(&(sample(0, 0))), &ldb, &beta, &prod(0, 0), &m, <ransb, <ransa);
416 }
417
418 return prod;
419 }
420
symVectProd(const Point & pt) const421 Point MatrixImplementation::symVectProd (const Point & pt) const
422 {
423 if (nbColumns_ != pt.getDimension() ) throw InvalidDimensionException(HERE) << "Invalid dimension in matrix/vector product";
424
425 Point prod(nbRows_);
426 // In this case, nbRows_ == nbColumns_
427 if (nbRows_ == 0) return prod;
428 char uplo('L');
429 int n(nbRows_);
430 int one(1);
431 double alpha(1.0);
432 double beta(0.0);
433 int luplo(1);
434 dsymv_(&uplo, &n, &alpha, const_cast<double*>(&((*this)[0])), &n, const_cast<double*>(&(pt[0])), &one, &beta, &prod[0], &one, &luplo);
435
436 return prod;
437 }
438
439 /* Gram matrix. If transpose == true, compute M^T.M, else M.M^T. */
computeGram(const Bool transposed) const440 MatrixImplementation MatrixImplementation::computeGram (const Bool transposed) const
441 {
442 if ((nbRows_ == 0) || (nbColumns_ == 0)) return MatrixImplementation(0, 0);
443 char uplo = 'L';
444 char trans = transposed ? 'T' : 'N';
445 int n = transposed ? nbColumns_ : nbRows_;
446 int k = transposed ? nbRows_ : nbColumns_;
447 double alpha = 1.0;
448 int lda = transposed ? k : n;
449 double beta = 0.0;
450 MatrixImplementation C(n, n);
451 int ldc = n;
452 int one = 1;
453 dsyrk_(&uplo, &trans, &n, &k, &alpha, const_cast<double*>(&((*this)[0])), &lda, &beta, &C[0], &ldc, &one, &one);
454
455 return C;
456 }
457
458 /* Multiplication by a Scalar */
operator *(const Scalar s) const459 MatrixImplementation MatrixImplementation::operator* (const Scalar s) const
460 {
461 if (s == 0.0) return MatrixImplementation(nbRows_, nbColumns_);
462 if ((nbRows_ == 0) || (nbColumns_ == 0)) return *this;
463 MatrixImplementation scalprod(*this);
464 double alpha(s);
465 int one(1);
466 int n_(nbRows_ * nbColumns_);
467 dscal_(&n_, &alpha, &scalprod[0], &one);
468
469 return scalprod;
470 }
471
472 /* In-place Multiplication by a Scalar */
operator *=(const Scalar s)473 MatrixImplementation & MatrixImplementation::operator*= (const Scalar s)
474 {
475 if ((nbRows_ == 0) || (nbColumns_ == 0)) return *this;
476 double alpha(s);
477 int one(1);
478 int n_(nbRows_ * nbColumns_);
479 dscal_(&n_, &alpha, &(*this)[0], &one);
480
481 return *this;
482 }
483
484 /* Division by a Scalar*/
operator /(const Scalar s) const485 MatrixImplementation MatrixImplementation::operator/ (const Scalar s) const
486 {
487 if (!(s < 0.0 || s > 0.0)) throw InvalidArgumentException(HERE);
488 if ((nbRows_ == 0) || (nbColumns_ == 0)) return *this;
489 MatrixImplementation scalprod(*this);
490 double alpha(1.0 / s);
491 int one(1);
492 int n_(nbRows_ * nbColumns_);
493 dscal_(&n_, &alpha, &scalprod[0], &one);
494
495 return scalprod;
496 }
497
498 /* In-place division by a Scalar */
operator /=(const Scalar s)499 MatrixImplementation & MatrixImplementation::operator/= (const Scalar s)
500 {
501 if (!(s < 0.0 || s > 0.0)) throw InvalidArgumentException(HERE);
502 if ((nbRows_ == 0) || (nbColumns_ == 0)) return *this;
503 double alpha(1.0 / s);
504 int one(1);
505 int n_(nbRows_ * nbColumns_);
506 dscal_(&n_, &alpha, &(*this)[0], &one);
507
508 return *this;
509 }
510
triangularProd(const MatrixImplementation & matrix,const char triangularSide,const char upperLower) const511 MatrixImplementation MatrixImplementation::triangularProd(const MatrixImplementation & matrix,
512 const char triangularSide,
513 const char upperLower) const
514 {
515 if (nbColumns_ != matrix.nbRows_) throw InvalidDimensionException(HERE) << "Invalid dimensions in matrix/matrix product left="
516 << nbRows_ << "x" << nbColumns_
517 << " right=" << matrix.nbRows_ << "x"
518 << matrix.nbColumns_;
519
520 MatrixImplementation mult(matrix);
521 if ((nbRows_ == 0) || (nbColumns_ == 0) || (matrix.nbColumns_ == 0)) return mult;
522 char side(triangularSide);
523 int lside = 1;
524 char uplo(upperLower);
525 int luplo = 1;
526 char trans('N');
527 int ltrans = 1;
528 char diag('N');
529 int ldiag = 1;
530 int m(nbRows_);
531 int n(matrix.nbColumns_);
532
533 Scalar alpha = 1.0;
534
535 // Lapack routine
536 dtrmm_(&side, &uplo, &trans, &diag, &m, &n, &alpha, const_cast<double*>(&((*this)[0])), &m, const_cast<double*>(&(mult[0])), &m, &lside, &luplo, <rans, &ldiag);
537 return mult;
538 }
539
540 /* Integer power, general matrix */
genPower(const UnsignedInteger n) const541 MatrixImplementation MatrixImplementation::genPower(const UnsignedInteger n) const
542 {
543 Bool first = true;
544 UnsignedInteger exponent = n;
545 MatrixImplementation y;
546 MatrixImplementation z(*this);
547 while (exponent > 0)
548 {
549 // t is the right bit of exponent
550 const UnsignedInteger t = exponent % 2;
551 // remove last bit from exponent
552 exponent /= 2;
553 // if right bit is 1
554 if (t != 0)
555 {
556 // if it is the rightmost bit equals to 1
557 if (first)
558 {
559 first = false;
560 y = z;
561 }
562 else y = y.genProd(z);
563 // if no more bit to consider
564 if (exponent == 0) return y;
565 }
566 // Square the contribution
567 z = z.genProd(z);
568 }
569 return y;
570 }
571
572 /* Integer power, symmetric matrix */
symPower(const UnsignedInteger n) const573 MatrixImplementation MatrixImplementation::symPower(const UnsignedInteger n) const
574 {
575 Bool first = true;
576 UnsignedInteger exponent = n;
577 MatrixImplementation y;
578 MatrixImplementation z(*this);
579 while (exponent > 0)
580 {
581 // t is the right bit of exponent
582 const UnsignedInteger t = exponent % 2;
583 // remove last bit from exponent
584 exponent /= 2;
585 // if right bit is 1
586 if (t != 0)
587 {
588 // if it is the rightmost bit equals to 1
589 if (first)
590 {
591 first = false;
592 y = z;
593 }
594 else y = y.symProd(z, 'L');
595 // if no more bit to consider
596 if (exponent == 0) return y;
597 }
598 // Square the contribution
599 z = z.symProd(z, 'L');
600 }
601 return y;
602 }
603
604 /* Empty returns true if there is no element in the MatrixImplementation */
isEmpty() const605 Bool MatrixImplementation::isEmpty() const
606 {
607 return ((nbRows_ == 0) || (nbColumns_ == 0) || (PersistentCollection<Scalar>::isEmpty()));
608 }
609
610 /* Returns true if triangular lower or upper */
isTriangular(Bool lower) const611 Bool MatrixImplementation::isTriangular(Bool lower) const
612 {
613 if ( nbRows_ == nbColumns_ )
614 {
615 for ( UnsignedInteger j = 1; j < nbColumns_; ++ j )
616 for ( UnsignedInteger i = 0; i < j; ++ i )
617 if ( std::abs( (*this)[lower ? convertPosition(i, j) : convertPosition(j, i)] ) > 0. )
618 return false;
619 return true;
620 }
621 else
622 return false;
623 }
624
625 /* Comparison operator */
operator ==(const MatrixImplementation & rhs) const626 Bool MatrixImplementation::operator == (const MatrixImplementation & rhs) const
627 {
628 const MatrixImplementation &lhs(*this);
629 Bool equality = true;
630
631 if (&lhs != &rhs) // Not the same object
632 {
633 const Collection<Scalar> & refLhs(static_cast<const Collection<Scalar> >(lhs));
634 const Collection<Scalar> & refRhs(static_cast<const Collection<Scalar> >(rhs));
635 equality = ( lhs.nbRows_ == rhs.nbRows_ && lhs.nbColumns_ == rhs.nbColumns_ && refLhs == refRhs);
636 }
637
638 return equality;
639 }
640
641
isSymmetric() const642 Bool MatrixImplementation::isSymmetric() const
643 {
644 const Scalar epsilon = ResourceMap::GetAsScalar("Matrix-SymmetryThreshold");
645 if ( nbRows_ == nbColumns_ )
646 {
647 for ( UnsignedInteger i = 1; i < nbRows_; ++ i )
648 for ( UnsignedInteger j = 0; j < i; ++ j )
649 if ( std::abs(this->operator[](convertPosition(i, j)) - operator[](convertPosition(j, i))) > epsilon)
650 return false;
651 return true;
652 }
653 else
654 return false;
655 }
656
657
658 /* Copy the lower triangle into the upper one, used by symmetric matrices */
symmetrize() const659 void MatrixImplementation::symmetrize() const
660 {
661 MatrixImplementation *refThis(const_cast<MatrixImplementation *>(this));
662 // Loop over the upper triangle
663 for (UnsignedInteger j = 0; j < nbColumns_; ++j)
664 for (UnsignedInteger i = 0; i < j; ++i)
665 refThis->operator[](convertPosition(i, j)) = operator[](convertPosition(j, i));
666 }
667
668 /* Fill the relevant triangle with zero */
triangularize(const Bool isLowerTriangular) const669 void MatrixImplementation::triangularize(const Bool isLowerTriangular) const
670 {
671 MatrixImplementation *refThis(const_cast<MatrixImplementation *>(this));
672 if (isLowerTriangular)
673 {
674 for (UnsignedInteger j = 0; j < nbColumns_; ++j)
675 for (UnsignedInteger i = 0; i < j; ++i)
676 refThis->operator[](convertPosition(i, j)) = 0.0;
677 }
678 else
679 {
680 for (UnsignedInteger j = 0; j < nbColumns_; ++j)
681 for (UnsignedInteger i = j + 1; i < nbRows_; ++i)
682 refThis->operator[](convertPosition(i, j)) = 0.0;
683 }
684 }
685
686 /* Check if the matrix values belong to (-1;1) */
hasUnitRange() const687 Bool MatrixImplementation::hasUnitRange() const
688 {
689 bool unitRange = true;
690 for (UnsignedInteger i = 0; i < nbRows_; ++i)
691 for (UnsignedInteger j = 0; j < nbColumns_; ++j)
692 {
693 if (std::abs(this->operator[](convertPosition(i, j))) > 1.0)
694 {
695 unitRange = false;
696 break;
697 }
698 }
699 return unitRange;
700 }
701
702 /* Set small elements to zero */
clean(const Scalar threshold) const703 MatrixImplementation MatrixImplementation::clean(const Scalar threshold) const
704 {
705 // Nothing to do for nonpositive threshold
706 if (threshold <= 0.0) return *this;
707 MatrixImplementation result(nbRows_, nbColumns_);
708 for (UnsignedInteger j = 0; j < nbColumns_; ++j)
709 for (UnsignedInteger i = 0; i < nbRows_; ++i)
710 {
711 const Scalar value = (*this)(i, j);
712 // Things are done this way to prevent spurious -0.0
713 if (std::abs(value) < 0.5 * threshold) result(i, j) = 0.0;
714 else result(i, j) = threshold * round(value / threshold);
715 }
716 return result;
717 }
718
719 /* Set small elements to zero */
cleanSym(const Scalar threshold) const720 MatrixImplementation MatrixImplementation::cleanSym(const Scalar threshold) const
721 {
722 symmetrize();
723 return clean(threshold);
724 }
725
726 /* Resolution of a linear system : rectangular matrix
727 * MX = b, M is an mxn matrix, b is an mxq matrix and
728 * X is an nxq matrix */
solveLinearSystemRect(const MatrixImplementation & b,const Bool keepIntact)729 MatrixImplementation MatrixImplementation::solveLinearSystemRect (const MatrixImplementation & b,
730 const Bool keepIntact)
731 {
732 if (nbRows_ != b.nbRows_) throw InvalidDimensionException(HERE) << "The right-hand side has row dimension=" << b.nbRows_ << ", expected " << nbRows_;
733 if (!(nbRows_ > 0 && nbColumns_ > 0 && b.nbColumns_ > 0)) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix or empty right-hand side";
734 int m(nbRows_);
735 int n(nbColumns_);
736 // B is an extended copy of b, it must be large enough to store the solution, see LAPACK documentation
737 int p(std::max(m, n));
738 int q(b.nbColumns_);
739 MatrixImplementation B(p, q);
740 for(UnsignedInteger j = 0; j < static_cast<UnsignedInteger>(q); ++j)
741 for (UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(m); ++i)
742 B(i, j) = b(i, j);
743 int nrhs(q);
744 int lwork(-1);
745 double lwork_d = -1.;
746 int info = -1;
747 std::vector<int> jpiv(n);
748 double rcond(ResourceMap::GetAsScalar("Matrix-DefaultSmallPivot"));
749 int rank = -1;
750
751 MatrixImplementation Q;
752 if (keepIntact) Q = MatrixImplementation(*this);
753 MatrixImplementation & A = keepIntact ? Q : *this;
754
755 dgelsy_(&m, &n, &nrhs, &A[0], &m, &B[0], &p, &jpiv[0], &rcond, &rank, &lwork_d, &lwork, &info);
756 lwork = static_cast<int>(lwork_d);
757 Point work(lwork);
758 dgelsy_(&m, &n, &nrhs, &A[0], &m, &B[0], &p, &jpiv[0], &rcond, &rank, &work[0], &lwork, &info);
759
760 MatrixImplementation result(n, q);
761 for(UnsignedInteger j = 0; j < static_cast<UnsignedInteger>(q); ++j)
762 for (UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(n); ++i)
763 result(i, j) = B(i, j);
764 return result;
765 }
766
767 /* Resolution of a linear system : rectangular matrix
768 * Mx = b, M is an mxn matrix, b is an m-dimensional
769 * vector and x is an n-dimensional vector */
solveLinearSystemRect(const Point & b,const Bool keepIntact)770 Point MatrixImplementation::solveLinearSystemRect (const Point & b,
771 const Bool keepIntact)
772 {
773 const UnsignedInteger m = b.getDimension();
774 if (nbRows_ != m) throw InvalidDimensionException(HERE) << "The right-hand side dimension is " << m << ", expected " << nbRows_;
775 if (nbRows_ == 0) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix";
776 // Solve the matrix linear system
777 // A MatrixImplementation is also a collection of Scalar, so it is automatically converted into a Point
778 return solveLinearSystemRect(MatrixImplementation(m, 1, b), keepIntact);
779 }
780
781
782 /* Resolution of a linear system : square matrix */
solveLinearSystemTri(const MatrixImplementation & b,const Bool keepIntact,const Bool lower,const Bool transposed)783 MatrixImplementation MatrixImplementation::solveLinearSystemTri (const MatrixImplementation & b,
784 const Bool keepIntact,
785 const Bool lower,
786 const Bool transposed)
787 {
788 if (nbRows_ != b.nbRows_ ) throw InvalidDimensionException(HERE) << "The right-hand side has row dimension=" << b.nbRows_ << ", expected " << nbRows_;
789 if (!(nbRows_ > 0 && nbColumns_ > 0 && b.nbColumns_ > 0)) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix or empty right-hand side";
790
791 // We must copy the right-hand side because it will be overwritten by the operation
792 MatrixImplementation B(b);
793 // side tells if we solve M.X = alpha.B or X.M = alpha.B
794 char side = 'L';
795 int lside = 1;
796 // M must be triangular. uplo tells if it is upper or lower triangular
797 char uplo = lower ? 'L' : 'U';
798 int luplo = 1;
799 // transa tells if M is transposed or not
800 char transa = transposed ? 'T' : 'N';
801 int ltransa = 1;
802 // diag tells if M is unit diagonal or not
803 char diag = 'N';
804 int ldiag = 1;
805 // the row dimension of M
806 int m = B.nbRows_;
807 // the column dimension of M
808 int n = B.nbColumns_;
809 // we solve the case alpha=1
810 double alpha = 1.0;
811 // leading dimension of M
812 int lda = nbRows_;
813 // leading dimension of B
814 int ldb = b.nbRows_;
815 if (keepIntact)
816 {
817 MatrixImplementation A(*this);
818 dtrsm_(&side, &uplo, &transa, &diag, &m, &n, &alpha, const_cast<double*>(&(A[0])), &lda, const_cast<double*>(&(B[0])), &ldb, &lside, &luplo, <ransa, &ldiag);
819 }
820 else dtrsm_(&side, &uplo, &transa, &diag, &m, &n, &alpha, const_cast<double*>(&((*this)[0])), &lda, const_cast<double*>(&(B[0])), &ldb, &lside, &luplo, <ransa, &ldiag);
821 return B;
822 }
823
824 /* Resolution of a linear system : square matrix */
solveLinearSystemTri(const Point & b,const Bool keepIntact,const Bool lower,const Bool transposed)825 Point MatrixImplementation::solveLinearSystemTri (const Point & b,
826 const Bool keepIntact,
827 const Bool lower,
828 const Bool transposed)
829 {
830 const UnsignedInteger m = b.getDimension();
831 if (nbRows_ != m) throw InvalidDimensionException(HERE) << "The right-hand side dimension is " << m << ", expected " << nbRows_;
832 if (nbRows_ == 0) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix";
833 // A MatrixImplementation is also a collection of Scalar, so it is automatically converted into a Point
834 return solveLinearSystemTri(MatrixImplementation(m, 1, b), keepIntact, lower, transposed);
835 }
836
837
838
839 /* Resolution of a linear system : square matrix */
solveLinearSystemSquare(const MatrixImplementation & b,const Bool keepIntact)840 MatrixImplementation MatrixImplementation::solveLinearSystemSquare (const MatrixImplementation & b,
841 const Bool keepIntact)
842 {
843 if (nbColumns_ != b.nbRows_ ) throw InvalidDimensionException(HERE) << "The right-hand side has row dimension=" << b.nbRows_ << ", expected " << nbRows_;
844 if (!(nbRows_ > 0 && nbColumns_ > 0 && b.nbColumns_ > 0)) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix or empty right-hand side";
845
846 // We must copy the right-hand side because it will be overwritten by the operation
847 MatrixImplementation B(b);
848 int m(nbRows_);
849 int nrhs(B.nbColumns_);
850 int info;
851 std::vector<int> ipiv(m);
852 if (keepIntact)
853 {
854 MatrixImplementation A(*this);
855 dgesv_(&m, &nrhs, &A[0], &m, &ipiv[0], &B[0], &m, &info);
856 }
857 else dgesv_(&m, &nrhs, &(*this)[0], &m, &ipiv[0], &B[0], &m, &info);
858 if (info != 0) throw NotDefinedException(HERE) << "Error: the matrix is singular.";
859 return B;
860 }
861
862 /* Resolution of a linear system : square matrix */
solveLinearSystemSquare(const Point & b,const Bool keepIntact)863 Point MatrixImplementation::solveLinearSystemSquare (const Point & b,
864 const Bool keepIntact)
865 {
866 const UnsignedInteger m = b.getDimension();
867 if (nbRows_ != m) throw InvalidDimensionException(HERE) << "The right-hand side dimension is " << m << ", expected " << nbRows_;
868 if (nbRows_ == 0) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix";
869 // A MatrixImplementation is also a collection of Scalar, so it is automatically converted into a Point
870 return solveLinearSystemRect(MatrixImplementation(m, 1, b), keepIntact);
871 }
872
873 /* Resolution of a linear system : symmetric matrix */
solveLinearSystemSym(const MatrixImplementation & b,const Bool keepIntact)874 MatrixImplementation MatrixImplementation::solveLinearSystemSym (const MatrixImplementation & b,
875 const Bool keepIntact)
876 {
877 if (nbColumns_ != b.nbRows_ ) throw InvalidDimensionException(HERE) << "The right-hand side has row dimension=" << b.nbRows_ << ", expected " << nbRows_;
878 if (!(nbRows_ > 0 && nbColumns_ > 0 && b.nbColumns_ > 0)) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix or empty right-hand side";
879
880 char uplo('L');
881 // We must copy the right-hand side as it will be overwritten by the operation
882 MatrixImplementation B(b);
883 int n(nbRows_);
884 int nrhs(B.nbColumns_);
885 int lwork(-1);
886 double lwork_d = -1.;
887 int info = -1;
888 std::vector<int> ipiv(n);
889 int luplo(1);
890
891 MatrixImplementation Q;
892 if (keepIntact) Q = MatrixImplementation(*this);
893 MatrixImplementation & A = keepIntact ? Q : *this;
894
895 dsysv_(&uplo, &n, &nrhs, &A[0], &n, &ipiv[0], &B[0], &n, &lwork_d, &lwork, &info, &luplo);
896 lwork = static_cast<int>(lwork_d);
897 Point work(lwork);
898 dsysv_(&uplo, &n, &nrhs, &A[0], &n, &ipiv[0], &B[0], &n, &work[0], &lwork, &info, &luplo);
899
900 if (info != 0) throw NotDefinedException(HERE) << "Error: the matrix is singular.";
901 return B;
902 }
903
904 /* Resolution of a linear system : symmetric matrix */
solveLinearSystemSym(const Point & b,const Bool keepIntact)905 Point MatrixImplementation::solveLinearSystemSym (const Point & b,
906 const Bool keepIntact)
907 {
908 const UnsignedInteger dimension = b.getDimension();
909 if (nbRows_ != dimension) throw InvalidDimensionException(HERE) << "The right-hand side dimension is " << dimension << ", expected " << nbRows_;
910 if (nbRows_ == 0) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix";
911 MatrixImplementation B(dimension, 1, b);
912 // A MatrixImplementation is also a collection of Scalar, so it is automatically converted into a Point
913 return solveLinearSystemSym(B, keepIntact);
914 }
915
916 /* Resolution of a linear system : covariance matrix */
solveLinearSystemCov(const MatrixImplementation & b,const Bool keepIntact)917 MatrixImplementation MatrixImplementation::solveLinearSystemCov (const MatrixImplementation & b,
918 const Bool keepIntact)
919 {
920 if (nbRows_ != b.nbRows_ ) throw InvalidDimensionException(HERE) << "The right-hand side has row dimension=" << b.nbRows_ << ", expected " << nbRows_;
921 if (!(nbRows_ > 0 && nbColumns_ > 0 && b.nbColumns_ > 0)) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix or empty right-hand side";
922
923 char uplo('L');
924 // We must copy the right-hand side as it will be overwritten by the operation
925 MatrixImplementation B(b);
926 int n(nbRows_);
927 int nrhs(B.nbColumns_);
928 int info;
929 int luplo(1);
930 if (keepIntact)
931 {
932 MatrixImplementation A(*this);
933 dposv_(&uplo, &n, &nrhs, &A[0], &n, &B[0], &n, &info, &luplo);
934 }
935 else
936 {
937 dposv_(&uplo, &n, &nrhs, &(*this)[0], &n, &B[0], &n, &info, &luplo);
938 }
939 if (info != 0) throw NotDefinedException(HERE) << "Error: the matrix is singular.";
940 return B;
941 }
942
943 /* Resolution of a linear system : symmetric matrix */
solveLinearSystemCov(const Point & b,const Bool keepIntact)944 Point MatrixImplementation::solveLinearSystemCov (const Point & b,
945 const Bool keepIntact)
946 {
947 const UnsignedInteger dimension = b.getDimension();
948 if (nbRows_ != dimension) throw InvalidDimensionException(HERE) << "The right-hand side dimension is " << dimension << ", expected " << nbRows_;
949 if (nbRows_ == 0) throw InvalidDimensionException(HERE) << "Cannot solve a linear system with empty matrix";
950 MatrixImplementation B(dimension, 1, b);
951 // A MatrixImplementation is also a collection of Scalar, so it is automatically converted into a Point
952 return solveLinearSystemCov(B, keepIntact);
953 }
954
955 /* Compute determinant */
computeLogAbsoluteDeterminant(Scalar & sign,const Bool keepIntact)956 Scalar MatrixImplementation::computeLogAbsoluteDeterminant (Scalar & sign,
957 const Bool keepIntact)
958 {
959 int n(nbRows_);
960 if (!(n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the determinant of an empty matrix";
961 Scalar logAbsoluteDeterminant = 0.0;
962 sign = 1.0;
963 if (n <= 2)
964 {
965 Scalar value = 0.0;
966 if (n == 1) value = (*this)(0, 0);
967 else value = (*this)(0, 0) * (*this)(1, 1) - (*this)(0, 1) * (*this)(1, 0);
968 if (value == 0.0)
969 {
970 sign = 0.0;
971 logAbsoluteDeterminant = SpecFunc::LowestScalar;
972 }
973 else
974 {
975 sign = (value > 0.0 ? 1.0 : -1.0);
976 logAbsoluteDeterminant = log(std::abs(value));
977 }
978 } // n <= 2
979 else
980 {
981 std::vector<int> ipiv (n);
982 int info = -1;
983
984 MatrixImplementation Q;
985 if (keepIntact) Q = MatrixImplementation(*this);
986 MatrixImplementation & A = keepIntact ? Q : *this;
987
988 // LU Factorization with LAPACK
989 dgetrf_(&n, &n, &A[0], &n, &ipiv[0], &info);
990 if (info > 0) return SpecFunc::LowestScalar;
991 // Determinant computation
992 for (UnsignedInteger i = 0; i < ipiv.size(); ++i)
993 {
994 const Scalar pivot = A[i * (ipiv.size() + 1)];
995 if (std::abs(pivot) == 0.0)
996 {
997 logAbsoluteDeterminant = SpecFunc::LowestScalar;
998 sign = 0.0;
999 }
1000 else
1001 {
1002 logAbsoluteDeterminant += log(std::abs(pivot));
1003 }
1004 if (pivot < 0.0) sign = -sign;
1005 if (ipiv[i] != int(i + 1)) sign = -sign;
1006 }
1007 } // n > 2
1008 return logAbsoluteDeterminant;
1009 }
1010
1011 /* Compute determinant */
computeDeterminant(const Bool keepIntact)1012 Scalar MatrixImplementation::computeDeterminant (const Bool keepIntact)
1013 {
1014 if (nbRows_ == 1) return (*this)(0, 0);
1015 if (nbRows_ == 2) return (*this)(0, 0) * (*this)(1, 1) - (*this)(0, 1) * (*this)(1, 0);
1016 Scalar sign = 0.0;
1017 const Scalar logAbsoluteDeterminant = computeLogAbsoluteDeterminant(sign, keepIntact);
1018 if (logAbsoluteDeterminant <= SpecFunc::LowestScalar) return 0.0;
1019 return sign * exp(logAbsoluteDeterminant);
1020 }
1021
1022 /* Compute determinant for a symmetric matrix */
computeLogAbsoluteDeterminantSym(Scalar & sign,const Bool keepIntact)1023 Scalar MatrixImplementation::computeLogAbsoluteDeterminantSym (Scalar & sign,
1024 const Bool keepIntact)
1025 {
1026 symmetrize();
1027 return computeLogAbsoluteDeterminant(sign, keepIntact);
1028 /* The implementation based on dsytrf does not uses the 2x2 diagonal blocks correctly
1029 int n(nbRows_);
1030 if (n == 0) throw InvalidDimensionException(HERE);
1031 std::vector<int> ipiv (n);
1032 char uplo('L');
1033 int info;
1034 Scalar determinant = 1.0;
1035 int lwork(-1);
1036 double lwork_d;
1037 int luplo(1);
1038
1039 // LU Factorization with LAPACK
1040 if (keepIntact)
1041 {
1042 MatrixImplementation A(*this);
1043 dsytrf_(&uplo, &n, &A[0], &n, &ipiv[0],&lwork_d, &lwork, &info, &luplo);
1044 lwork = static_cast<int>(lwork_d);
1045 Point work(lwork);
1046 dsytrf_(&uplo, &n, &A[0], &n, &ipiv[0],&work[0], &lwork, &info, &luplo);
1047 // Determinant computation
1048 for (UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(n); ++i)
1049 {
1050 determinant *= A[i * (ipiv.size() + 1)];
1051 if (ipiv[i] != int(i + 1)) determinant = -determinant;
1052 }
1053 }
1054 else
1055 {
1056 dsytrf_(&uplo, &n, &(*this)[0], &n, &ipiv[0],&lwork_d, &lwork, &info, &luplo);
1057 lwork = static_cast<int>(lwork_d);
1058 Point work(lwork);
1059 dsytrf_(&uplo, &n, &(*this)[0], &n, &ipiv[0],&work[0], &lwork, &info, &luplo);
1060 // Determinant computation
1061 for (UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(n); ++i)
1062 {
1063 determinant *= (*this)[i * (ipiv.size() + 1)];
1064 if (ipiv[i] != int(i + 1)) determinant = -determinant;
1065 }
1066 }
1067
1068 return determinant; */
1069 }
1070
1071 /* Compute determinant for a symmetric matrix */
computeDeterminantSym(const Bool keepIntact)1072 Scalar MatrixImplementation::computeDeterminantSym (const Bool keepIntact)
1073 {
1074 if (nbRows_ == 1) return (*this)(0, 0);
1075 if (nbRows_ == 2) return (*this)(0, 0) * (*this)(1, 1) - (*this)(1, 0) * (*this)(1, 0);
1076 Scalar sign = 0.0;
1077 const Scalar logAbsoluteDeterminant = computeLogAbsoluteDeterminant(sign, keepIntact);
1078 if (logAbsoluteDeterminant <= SpecFunc::LowestScalar) return 0.0;
1079 return sign * exp(logAbsoluteDeterminant);
1080 }
1081
1082 /* Compute trace */
computeTrace() const1083 Scalar MatrixImplementation::computeTrace() const
1084 {
1085 Scalar trace = 0.0;
1086 for (UnsignedInteger i = 0; i < nbRows_; ++i) trace += (*this)(i, i);
1087 return trace;
1088 }
1089
1090 /* Compute the eigenvalues of a square matrix */
computeEigenValuesSquare(const Bool keepIntact)1091 MatrixImplementation::ComplexCollection MatrixImplementation::computeEigenValuesSquare (const Bool keepIntact)
1092 {
1093 int n(nbRows_);
1094 if (!(n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the eigenvalues of an empty matrix";
1095 char jobvl('N');
1096 char jobvr('N');
1097 Point wr(n, 0.0);
1098 Point wi(n, 0.0);
1099 double vl = 0.;
1100 double vr = 0.;
1101 int ldvl(1);
1102 int ldvr(1);
1103 int lwork(-1);
1104 double lwork_d = -1.;
1105 int info = -1;
1106 int ljobvl(1);
1107 int ljobvr(1);
1108
1109 MatrixImplementation Q;
1110 if (keepIntact) Q = MatrixImplementation(*this);
1111 MatrixImplementation & A = keepIntact ? Q : *this;
1112
1113 dgeev_(&jobvl, &jobvr, &n, &A[0], &n, &wr[0], &wi[0], &vl, &ldvl, &vr, &ldvr, &lwork_d, &lwork, &info, &ljobvl, &ljobvr);
1114 lwork = static_cast<int>(lwork_d);
1115 Point work(lwork);
1116 dgeev_(&jobvl, &jobvr, &n, &A[0], &n, &wr[0], &wi[0], &vl, &ldvl, &vr, &ldvr, &work[0], &lwork, &info, &ljobvl, &ljobvr);
1117
1118 if (info != 0) throw InternalException(HERE) << "Error: the QR algorithm failed to converge.";
1119 ComplexCollection eigenValues(n);
1120 for (UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(n); ++i) eigenValues[i] = Complex(wr[i], wi[i]);
1121 return eigenValues;
1122 }
1123
computeEVSquare(ComplexMatrixImplementation & v,const Bool keepIntact)1124 MatrixImplementation::ComplexCollection MatrixImplementation::computeEVSquare (ComplexMatrixImplementation & v,
1125 const Bool keepIntact)
1126 {
1127 int n(nbRows_);
1128 if (!(n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the eigenvalues of an empty matrix";
1129 char jobvl('N');
1130 char jobvr('V');
1131 Point wr(n, 0.0);
1132 Point wi(n, 0.0);
1133 double vl;
1134 MatrixImplementation vr(n, n);
1135 int ldvl(1);
1136 int ldvr(n);
1137 int lwork(-1);
1138 double lwork_d;
1139 int info;
1140 int ljobvl(1);
1141 int ljobvr(1);
1142
1143 MatrixImplementation Q;
1144 if (keepIntact) Q = MatrixImplementation(*this);
1145 MatrixImplementation & A = keepIntact ? Q : *this;
1146
1147 dgeev_(&jobvl, &jobvr, &n, &A[0], &n, &wr[0], &wi[0], &vl, &ldvl, &vr[0], &ldvr, &lwork_d, &lwork, &info, &ljobvl, &ljobvr);
1148 lwork = static_cast<int>(lwork_d);
1149 Point work(lwork);
1150 dgeev_(&jobvl, &jobvr, &n, &A[0], &n, &wr[0], &wi[0], &vl, &ldvl, &vr[0], &ldvr, &work[0], &lwork, &info, &ljobvl, &ljobvr);
1151
1152 // Cast the eigenvalues into OpenTURNS data structures
1153 ComplexCollection eigenValues(n);
1154 for (UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(n); ++i)
1155 {
1156 eigenValues[i] = Complex(wr[i], wi[i]);
1157 }
1158 if (info != 0) throw InternalException(HERE) << "Error: the QR algorithm failed to converge.";
1159 // Cast the eigenvectors into OpenTURNS data structures
1160 v = ComplexMatrixImplementation(n, n);
1161 UnsignedInteger j = 0;
1162 while (j < static_cast<UnsignedInteger>(n))
1163 {
1164 // Real eigenvalue
1165 if (wi[j] == 0.0)
1166 {
1167 for (UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(n); ++i) v(i, j) = vr(i, j);
1168 ++j;
1169 }
1170 // Complex conjugate pair of eigenvalues
1171 else
1172 {
1173 for (UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(n); ++i)
1174 {
1175 v(i, j) = Complex(vr(i, j), vr(i, j + 1));
1176 v(i, j + 1) = Complex(vr(i, j), -vr(i, j + 1));
1177 }
1178 j += 2;
1179 }
1180 } // Conversion of eigenvectors
1181 return eigenValues;
1182 }
1183
1184 /* Compute the eigenvalues of a symmetric matrix */
computeEigenValuesSym(const Bool keepIntact)1185 Point MatrixImplementation::computeEigenValuesSym (const Bool keepIntact)
1186 {
1187 int n(nbRows_);
1188 if (!(n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the eigenvalues of an empty matrix";
1189 char jobz('N');
1190 char uplo('L');
1191 Point w(n, 0.0);
1192 int lwork(-1);
1193 double lwork_d;
1194 int info;
1195 int ljobz(1);
1196 int luplo(1);
1197
1198 MatrixImplementation Q;
1199 if (keepIntact) Q = MatrixImplementation(*this);
1200 MatrixImplementation & A = keepIntact ? Q : *this;
1201
1202 dsyev_(&jobz, &uplo, &n, &A[0], &n, &w[0], &lwork_d, &lwork, &info, &ljobz, &luplo);
1203 lwork = static_cast<int>(lwork_d);
1204 Point work(lwork);
1205 dsyev_(&jobz, &uplo, &n, &A[0], &n, &w[0], &work[0], &lwork, &info, &ljobz, &luplo);
1206
1207 if (info != 0) throw InternalException(HERE) << "Error: the QR algorithm failed to converge.";
1208 return w;
1209 }
1210
computeEVSym(MatrixImplementation & v,const Bool keepIntact)1211 Point MatrixImplementation::computeEVSym (MatrixImplementation & v,
1212 const Bool keepIntact)
1213 {
1214 int n(nbRows_);
1215 if (!(n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the eigenvalues of an empty matrix";
1216 char jobz('V');
1217 char uplo('L');
1218 Point w(n, 0.0);
1219 int lwork(-1);
1220 double lwork_d;
1221 int info;
1222 int ljobz(1);
1223 int luplo(1);
1224
1225 MatrixImplementation Q;
1226 if (keepIntact) Q = MatrixImplementation(*this);
1227 MatrixImplementation & A = keepIntact ? Q : *this;
1228
1229 dsyev_(&jobz, &uplo, &n, &A[0], &n, &w[0], &lwork_d, &lwork, &info, &ljobz, &luplo);
1230 lwork = static_cast<int>(lwork_d);
1231 Point work(lwork);
1232 dsyev_(&jobz, &uplo, &n, &A[0], &n, &w[0], &work[0], &lwork, &info, &ljobz, &luplo);
1233 v = A;
1234
1235 if (info != 0) throw InternalException(HERE) << "Error: the QR algorithm failed to converge.";
1236 return w;
1237 }
1238
1239 /* Compute the largest eigenvalue module using power iterations, square matrix */
computeLargestEigenValueModuleSquare(Scalar & maximumModule,const UnsignedInteger maximumIterations,const Scalar epsilon) const1240 Bool MatrixImplementation::computeLargestEigenValueModuleSquare(Scalar & maximumModule,
1241 const UnsignedInteger maximumIterations,
1242 const Scalar epsilon) const
1243 {
1244 const UnsignedInteger dimension = getNbRows();
1245 Point currentEigenVector(dimension, 1.0);
1246 Point nextEigenVector(genVectProd(currentEigenVector));
1247 Scalar nextEigenValue = nextEigenVector.norm();
1248 maximumModule = nextEigenValue / std::sqrt(1.0 * dimension);
1249 Bool found = false;
1250 Scalar precision = 0.0;
1251 for (UnsignedInteger iteration = 0; iteration < maximumIterations && !found; ++iteration)
1252 {
1253 LOGDEBUG(OSS() << "(" << iteration << ") maximum module=" << maximumModule);
1254 currentEigenVector = nextEigenVector / nextEigenValue;
1255 nextEigenVector = genVectProd(currentEigenVector);
1256 nextEigenValue = nextEigenVector.norm();
1257 precision = std::abs(nextEigenValue - maximumModule);
1258 found = precision <= epsilon * nextEigenValue;
1259 LOGDEBUG(OSS() << "(" << iteration << ") precison=" << precision << ", relative precision=" << precision / nextEigenValue << ", found=" << found);
1260 maximumModule = nextEigenValue;
1261 }
1262 return found;
1263 }
1264
1265 /* Compute the largest eigenvalue module using power iterations, symmetric matrix */
computeLargestEigenValueModuleSym(Scalar & maximumModule,const UnsignedInteger maximumIterations,const Scalar epsilon) const1266 Bool MatrixImplementation::computeLargestEigenValueModuleSym(Scalar & maximumModule,
1267 const UnsignedInteger maximumIterations,
1268 const Scalar epsilon) const
1269 {
1270 const UnsignedInteger dimension = getNbRows();
1271 Point currentEigenVector(dimension, 1.0);
1272 Point nextEigenVector(symVectProd(currentEigenVector));
1273 Scalar nextEigenValue = nextEigenVector.norm();
1274 maximumModule = nextEigenValue / std::sqrt(1.0 * dimension);
1275 Bool found = false;
1276 Scalar precision = 0.0;
1277 for (UnsignedInteger iteration = 0; iteration < maximumIterations && !found; ++iteration)
1278 {
1279 LOGDEBUG(OSS() << "(" << iteration << ") maximum module=" << maximumModule);
1280 currentEigenVector = nextEigenVector / nextEigenValue;
1281 nextEigenVector = symVectProd(currentEigenVector);
1282 nextEigenValue = nextEigenVector.norm();
1283 precision = std::abs(nextEigenValue - maximumModule);
1284 found = precision <= epsilon * nextEigenValue;
1285 LOGDEBUG(OSS() << "(" << iteration << ") precison=" << precision << ", relative precision=" << precision / nextEigenValue << ", found=" << found);
1286 maximumModule = nextEigenValue;
1287 }
1288 return found;
1289 }
1290
1291 /* Compute the singular values of a matrix */
computeSingularValues(const Bool keepIntact)1292 Point MatrixImplementation::computeSingularValues(const Bool keepIntact)
1293 {
1294 int m = nbRows_;
1295 int n = nbColumns_;
1296 if (!(m > 0 && n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the singular values of an empty matrix";
1297
1298 // check for nans, cf https://github.com/Reference-LAPACK/lapack/issues/469
1299 for (UnsignedInteger j = 0; j < getNbColumns(); ++ j)
1300 for (UnsignedInteger i = 0; i < getNbRows(); ++ i)
1301 if (!SpecFunc::IsNormal(operator()(i, j)))
1302 throw InvalidArgumentException(HERE) << "Cannot compute singular values due to nan/inf values";
1303
1304 char jobz = 'N';
1305 Point S(std::min(m, n), 0.0);
1306 Point work(1, 0.0);
1307 MatrixImplementation u(1, 1);
1308 int ldu = 1;
1309 int ldvt = 1;
1310 MatrixImplementation vT(1, 1);
1311 int lwork = -1;
1312 std::vector<int> iwork(8 * std::min(m, n));
1313 int info = 0;
1314 int ljobz = 1;
1315
1316 MatrixImplementation Q;
1317 if (keepIntact) Q = MatrixImplementation(*this);
1318 MatrixImplementation & A = keepIntact ? Q : *this;
1319
1320 // First call to compute the optimal work size
1321 dgesdd_(&jobz, &m, &n, &A[0], &m, &S[0], &u[0], &ldu, &vT[0], &ldvt, &work[0], &lwork, &iwork[0], &info, &ljobz);
1322 lwork = static_cast<int>(work[0]);
1323 work = Point(lwork, 0.0);
1324 // Second call to compute the SVD
1325 dgesdd_(&jobz, &m, &n, &A[0], &m, &S[0], &u[0], &ldu, &vT[0], &ldvt, &work[0], &lwork, &iwork[0], &info, &ljobz);
1326
1327 if (info != 0) throw InternalException(HERE) << "Error: the updating process failed.";
1328 return S;
1329 }
1330
1331 /* Compute the singular values and singular decomposition of a matrix */
computeSVD(MatrixImplementation & u,MatrixImplementation & vT,const Bool fullSVD,const Bool keepIntact)1332 Point MatrixImplementation::computeSVD(MatrixImplementation & u,
1333 MatrixImplementation & vT,
1334 const Bool fullSVD,
1335 const Bool keepIntact)
1336 {
1337 int m(nbRows_);
1338 int n(nbColumns_);
1339 if (!(m > 0 && n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the singular values decomposition of an empty matrix";
1340 char jobz( fullSVD ? 'A' : 'S');
1341 int ldu(m);
1342 u = MatrixImplementation(m, ( fullSVD ? m : std::min(m, n)));
1343 int ldvt = (fullSVD ? n : std::min(m, n));
1344 vT = MatrixImplementation(( fullSVD ? n : std::min(m, n)), n);
1345 Point S(std::min(m, n), 0.0);
1346 Point work(1, 0.0);
1347 int lwork(-1);
1348 std::vector<int> iwork(8 * std::min(m, n));
1349 int info(0);
1350 int ljobz(1);
1351
1352 MatrixImplementation Q;
1353 if (keepIntact) Q = MatrixImplementation(*this);
1354 MatrixImplementation & A = keepIntact ? Q : *this;
1355
1356 // First call to compute the optimal work size
1357 dgesdd_(&jobz, &m, &n, &A[0], &m, &S[0], &u[0], &ldu, &vT[0], &ldvt, &work[0], &lwork, &iwork[0], &info, &ljobz);
1358 lwork = static_cast<int>(work[0]);
1359 work = Point(lwork, 0.0);
1360 // Second call to compute the SVD
1361 dgesdd_(&jobz, &m, &n, &A[0], &m, &S[0], &u[0], &ldu, &vT[0], &ldvt, &work[0], &lwork, &iwork[0], &info, &ljobz);
1362
1363 if (info != 0) throw InternalException(HERE) << "Error: LAPACK trouble in computing SVD decomposition, return code is " << info;
1364 return S;
1365 }
1366
1367
1368 /* Check if the matrix is SPD */
isPositiveDefinite() const1369 Bool MatrixImplementation::isPositiveDefinite() const
1370 {
1371 int info;
1372 int n(nbRows_);
1373 if (!(n > 0)) throw InvalidDimensionException(HERE) << "Cannot check the definite positiveness of an empty matrix";
1374 char uplo('L');
1375 int luplo(1);
1376 MatrixImplementation A(*this);
1377 dpotrf_(&uplo, &n, &A[0], &n, &info, &luplo);
1378 return (info == 0) ;
1379 }
1380
triangularVectProd(const ScalarCollection & pt,const char side,const Bool transpose) const1381 MatrixImplementation::ScalarCollection MatrixImplementation::triangularVectProd(const ScalarCollection & pt,
1382 const char side,
1383 const Bool transpose) const
1384 {
1385 char uplo = side;
1386 int luplo = 1;
1387
1388 // trans tells if the matrix is transposed or not
1389 char trans = transpose ? 'T' : 'N';
1390 int ltrans = 1;
1391
1392 // diag tells if M is unit diagonal or not
1393 char diag = 'N';
1394 int ldiag = 1;
1395
1396 // the dimension of the matrix
1397 int n = nbRows_;
1398
1399 // leading dimension of M
1400 int lda = nbRows_;
1401 int one = 1;
1402
1403 ScalarCollection x(nbRows_);
1404 for (UnsignedInteger i = 0; i < pt.getSize(); ++i) x[i] = pt[i];
1405
1406 dtrmv_(&uplo, &trans, &diag, &n, const_cast<double*>(&((*this)[0])), &lda, const_cast<double*>(&(x[0])), &one, &luplo, <rans, &ldiag);
1407 return x;
1408 }
1409
triangularVectProd(const Point & pt,const char side,const Bool transpose) const1410 MatrixImplementation::ScalarCollection MatrixImplementation::triangularVectProd(const Point & pt,
1411 const char side,
1412 const Bool transpose) const
1413 {
1414 return triangularVectProd(pt.getCollection(), side, transpose);
1415 }
1416
1417 /* Build the Cholesky factorization of the matrix */
computeCholesky(const Bool keepIntact)1418 MatrixImplementation MatrixImplementation::computeCholesky(const Bool keepIntact)
1419 {
1420 int n = nbRows_;
1421 if (!(n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the Cholesky decomposition of an empty matrix";
1422 int info = 0;
1423 char uplo = 'L';
1424 int luplo = 1;
1425
1426 MatrixImplementation L;
1427 if (keepIntact) L = MatrixImplementation(*this);
1428 MatrixImplementation & A = keepIntact ? L : *this;
1429
1430 dpotrf_(&uplo, &n, &A[0], &n, &info, &luplo);
1431 if (info != 0) throw NotSymmetricDefinitePositiveException(HERE) << "Error: the matrix is not definite positive.";
1432 for (UnsignedInteger j = 0; j < (UnsignedInteger)(n); ++ j)
1433 for (UnsignedInteger i = 0; i < j; ++ i)
1434 A(i, j) = 0.0;
1435 A.setName("");
1436 return A;
1437 }
1438
1439 /* 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 */
CholeskyUpdate(MatrixImplementation & cholesky,const Point & vector)1440 void MatrixImplementation::CholeskyUpdate(MatrixImplementation & cholesky,
1441 const Point & vector)
1442 {
1443 UnsignedInteger dimension = cholesky.nbRows_;
1444 if (dimension != cholesky.nbColumns_) throw InvalidDimensionException(HERE) << "Cannot update a non-square Cholesky factor";
1445 if (dimension != vector.getDimension()) throw InvalidDimensionException(HERE) << "Incompatible Cholesky factor dimension and vector dimension";
1446 // Working copy of vector
1447 Point work(vector);
1448 int size = dimension;
1449 int one = 1;
1450 UnsignedInteger shift = 0;
1451 // Parameters of the Givens rotation
1452 double cosI = 0.;
1453 double sinI = 0.;
1454 for (UnsignedInteger i = 0; i < dimension - 1; ++i)
1455 {
1456 // Generate Givens rotation
1457 drotg_(&cholesky[shift], &work[i], &cosI, &sinI);
1458 // Flip rotation if negative diagonal entry
1459 if (cholesky[shift] < 0.0)
1460 {
1461 cholesky[shift] = -cholesky[shift];
1462 cosI = -cosI;
1463 sinI = -sinI;
1464 }
1465 --size;
1466 // Perform the rotation
1467 drot_(&size, &cholesky[shift + 1], &one, &work[i + 1], &one, &cosI, &sinI);
1468 shift += dimension + 1;
1469 }
1470 // Last rotation
1471 // Generate Givens rotation
1472 drotg_(&cholesky[shift], &work[dimension - 1], &cosI, &sinI);
1473 // Flip rotation if negative diagonal entry
1474 if (cholesky[shift] < 0.0) cholesky[shift] = -cholesky[shift];
1475 }
1476
1477 /* 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 */
CholeskyDowndate(MatrixImplementation & cholesky,const Point & vector)1478 void MatrixImplementation::CholeskyDowndate(MatrixImplementation & cholesky,
1479 const Point & vector)
1480 {
1481 UnsignedInteger dimension = cholesky.nbRows_;
1482 if (dimension != cholesky.nbColumns_) throw InvalidDimensionException(HERE) << "Cannot update a non-square Cholesky factor";
1483 if (dimension != vector.getDimension()) throw InvalidDimensionException(HERE) << "Incompatible Cholesky factor dimension and vector dimension";
1484 // Working copy of vector
1485 Point work(cholesky.solveLinearSystemTri(vector));
1486 // Parameters of the Givens rotation
1487 Scalar qs = sqrt(1.0 - work.normSquare());
1488 Point cosI(dimension);
1489 Point sinI(dimension);
1490 for (SignedInteger i = dimension - 1; i >= 0; --i)
1491 {
1492 // Generate Givens rotation
1493 drotg_(&qs, &work[i], &cosI[i], &sinI[i]);
1494 // Flip rotation if negative diagonal entry
1495 if (qs < 0.0)
1496 {
1497 qs = -qs;
1498 cosI[i] = -cosI[i];
1499 sinI[i] = -sinI[i];
1500 }
1501 }
1502 work = Point(dimension, 0.0);
1503 int size(0);
1504 int one(1);
1505 UnsignedInteger shift = (dimension - 1) * (dimension + 1);
1506 for (SignedInteger i = dimension - 1; i >= 0; --i)
1507 {
1508 ++size;
1509 // Perform the rotation
1510 drot_(&size, &work[i], &one, &cholesky[shift], &one, &cosI[i], &sinI[i]);
1511 if (cholesky[shift] < 0.0)
1512 {
1513 qs = -1.0;
1514 dscal_(&size, &qs, &cholesky[shift], &one);
1515 }
1516 shift -= dimension + 1;
1517 }
1518 }
1519
1520 /* Build the QR factorization of the mxn matrix A
1521 Case m >= n:
1522 [R]
1523 A = [Q1 Q2][0]
1524 with Q1 a mxn matrix, Q2 a (m-n)xm matrix, R a nxn upper triangular matrix.
1525 [R]
1526 If fullQR == true, the matrices [Q1 Q2] and [0] are returned
1527 Else it is the matrices [Q1] and [R] that are returned.
1528 Case m < n:
1529 A = [Q][R1 R2]
1530 with Q a mxm matrix, R1 a mxm upper triangular matrix, R2 a mx(m-n) matrix.
1531 If fullQR == true or false, the matrices [Q] and [R1 R2] are returned
1532 */
computeQR(MatrixImplementation & R,const Bool fullQR,const Bool keepIntact)1533 MatrixImplementation MatrixImplementation::computeQR(MatrixImplementation & R,
1534 const Bool fullQR,
1535 const Bool keepIntact)
1536 {
1537 int m = nbRows_;
1538 int n = nbColumns_;
1539 int lda = nbRows_;
1540
1541 if (!(m > 0 && n > 0)) throw InvalidDimensionException(HERE) << "Cannot compute the QR decomposition of an empty matrix";
1542 int k = std::min(m, n);
1543 Point tau(k);
1544 int lwork = -1;
1545 int info = -1;
1546 double lwork_d = -1.;
1547 MatrixImplementation A;
1548 if (keepIntact) A = MatrixImplementation(*this);
1549 MatrixImplementation & Q = keepIntact ? A : *this;
1550
1551 // First call to dgeqrf to get the optimal size for the working space
1552 dgeqrf_(&m, &n, &Q[0], &lda, &tau[0], &lwork_d, &lwork, &info);
1553 if (info != 0) throw InternalException(HERE) << "Lapack DGEQRF: error code=" << info;
1554 // Here is the optimal size of the working space
1555 lwork = static_cast<int>(lwork_d);
1556 Point work(lwork);
1557 // Second call to compute the decomposition
1558 dgeqrf_(&m, &n, &Q[0], &lda, &tau[0], &work[0], &lwork, &info);
1559 if (info != 0) throw InternalException(HERE) << "Lapack DGEQRF: error code=" << info;
1560
1561 // Rebuild R
1562 int p(fullQR ? m : k);
1563 R = MatrixImplementation(p, n);
1564 for ( UnsignedInteger i = 0; i < static_cast<UnsignedInteger>(k) ; ++ i )
1565 for ( UnsignedInteger j = i; j < static_cast<UnsignedInteger>(n); ++ j )
1566 R(i, j) = Q(i, j);
1567
1568 // Rebuild Q
1569 // It is done using the product of the reflectors computed by dgeqrf
1570 // First call to dorgqr to get the optimal size for the working space
1571 lwork = -1;
1572 dorgqr_(&m, &p, &k, &Q[0], &lda, &tau[0], &lwork_d, &lwork, &info);
1573 if (info != 0) throw InternalException(HERE) << "Lapack DORGQR: error code=" << info;
1574 // Here is the optimal size of the working space
1575 lwork = static_cast<int>(lwork_d);
1576 work = Point(lwork);
1577 // Second call to compute the product of reflectors
1578 if (fullQR && (m > n))
1579 {
1580 // Here we must copy Q into a larger matrix to get the desired mxm Q factor before the call to dorgqr
1581 Q.resize(m * m);
1582 Q.nbRows_ = m;
1583 Q.nbColumns_ = m;
1584 }
1585 dorgqr_(&m, &p, &k, &Q[0], &lda, &tau[0], &work[0], &lwork, &info);
1586 if (m < n)
1587 {
1588 // Here we must copy Q into a smaller matrix to get the desired mxm Q factor after the call to dorgqr
1589 Q.resize(m * m);
1590 Q.nbRows_ = m;
1591 Q.nbColumns_ = m;
1592 }
1593 if (info != 0) throw InternalException(HERE) << "Lapack DORGQR: error code=" << info;
1594 Q.setName("");
1595 return Q;
1596 }
1597
1598 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const1599 void MatrixImplementation::save(Advocate & adv) const
1600 {
1601 PersistentCollection<Scalar>::save(adv);
1602 adv.saveAttribute( "nbRows_", nbRows_);
1603 adv.saveAttribute( "nbColumns_", nbColumns_);
1604 }
1605
1606 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)1607 void MatrixImplementation::load(Advocate & adv)
1608 {
1609 PersistentCollection<Scalar>::load(adv);
1610
1611 adv.loadAttribute( "nbRows_", nbRows_);
1612 adv.loadAttribute( "nbColumns_", nbColumns_);
1613 }
1614
stride(const UnsignedInteger dim) const1615 UnsignedInteger MatrixImplementation::stride(const UnsignedInteger dim) const
1616 {
1617 UnsignedInteger stride = elementSize();
1618 if (dim > 0)
1619 stride *= nbRows_;
1620 return stride;
1621 }
1622
1623
1624 END_NAMESPACE_OPENTURNS
1625