1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
5 Copyright (C) 2003, 2004, 2005, 2006 StatPro Italia srl
6 Copyright (C) 2003, 2004 Ferdinando Ametrano
7 Copyright (C) 2015 Michael von den Driesch
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22/*! \file matrix.hpp
23 \brief matrix used in linear algebra.
24*/
25
26#ifndef quantlib_matrix_hpp
27#define quantlib_matrix_hpp
28
29#include <ql/math/array.hpp>
30#include <ql/utilities/steppingiterator.hpp>
31#include <initializer_list>
32#include <iterator>
33
34namespace QuantLib {
35
36 //! %Matrix used in linear algebra.
37 /*! This class implements the concept of Matrix as used in linear
38 algebra. As such, it is <b>not</b> meant to be used as a
39 container.
40 */
41 class Matrix {
42 public:
43 //! \name Constructors, destructor, and assignment
44 //@{
45 //! creates a null matrix
46 Matrix();
47 //! creates a matrix with the given dimensions
48 Matrix(Size rows, Size columns);
49 //! creates the matrix and fills it with <tt>value</tt>
50 Matrix(Size rows, Size columns, Real value);
51 //! creates the matrix and fills it with data from a range.
52 /*! \warning if the range defined by [begin, end) is larger
53 than the size of the matrix, a memory access violation
54 might occur. It is up to the user to avoid this.
55 */
56 template <class Iterator>
57 Matrix(Size rows, Size columns, Iterator begin, Iterator end);
58 Matrix(const Matrix&);
59 Matrix(Matrix&&) noexcept;
60 Matrix(std::initializer_list<std::initializer_list<Real>>);
61 ~Matrix() = default;
62
63 Matrix& operator=(const Matrix&);
64 Matrix& operator=(Matrix&&) noexcept;
65
66 bool operator==(const Matrix&) const;
67 bool operator!=(const Matrix&) const;
68 //@}
69
70 //! \name Algebraic operators
71 /*! \pre all matrices involved in an algebraic expression must have
72 the same size.
73 */
74 //@{
75 const Matrix& operator+=(const Matrix&);
76 const Matrix& operator-=(const Matrix&);
77 const Matrix& operator*=(Real);
78 const Matrix& operator/=(Real);
79 //@}
80
81 typedef Real* iterator;
82 typedef const Real* const_iterator;
83 typedef std::reverse_iterator<iterator> reverse_iterator;
84 typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
85 typedef Real* row_iterator;
86 typedef const Real* const_row_iterator;
87 typedef std::reverse_iterator<row_iterator> reverse_row_iterator;
88 typedef std::reverse_iterator<const_row_iterator>
89 const_reverse_row_iterator;
90 typedef step_iterator<iterator> column_iterator;
91 typedef step_iterator<const_iterator> const_column_iterator;
92 typedef std::reverse_iterator<column_iterator>
93 reverse_column_iterator;
94 typedef std::reverse_iterator<const_column_iterator>
95 const_reverse_column_iterator;
96 //! \name Iterator access
97 //@{
98 const_iterator begin() const;
99 iterator begin();
100 const_iterator end() const;
101 iterator end();
102 const_reverse_iterator rbegin() const;
103 reverse_iterator rbegin();
104 const_reverse_iterator rend() const;
105 reverse_iterator rend();
106 const_row_iterator row_begin(Size i) const;
107 row_iterator row_begin(Size i);
108 const_row_iterator row_end(Size i) const;
109 row_iterator row_end(Size i);
110 const_reverse_row_iterator row_rbegin(Size i) const;
111 reverse_row_iterator row_rbegin(Size i);
112 const_reverse_row_iterator row_rend(Size i) const;
113 reverse_row_iterator row_rend(Size i);
114 const_column_iterator column_begin(Size i) const;
115 column_iterator column_begin(Size i);
116 const_column_iterator column_end(Size i) const;
117 column_iterator column_end(Size i);
118 const_reverse_column_iterator column_rbegin(Size i) const;
119 reverse_column_iterator column_rbegin(Size i);
120 const_reverse_column_iterator column_rend(Size i) const;
121 reverse_column_iterator column_rend(Size i);
122 //@}
123
124 //! \name Element access
125 //@{
126 const_row_iterator operator[](Size) const;
127 const_row_iterator at(Size) const;
128 row_iterator operator[](Size);
129 row_iterator at(Size);
130 Array diagonal() const;
131 Real& operator()(Size i, Size j) const;
132 //@}
133
134 //! \name Inspectors
135 //@{
136 Size rows() const;
137 Size columns() const;
138 bool empty() const;
139 Size size1() const;
140 Size size2() const;
141 //@}
142
143 //! \name Utilities
144 //@{
145 void swap(Matrix&) noexcept;
146 //@}
147 private:
148 std::unique_ptr<Real[]> data_;
149 Size rows_ = 0, columns_ = 0;
150 };
151
152 // algebraic operators
153
154 /*! \relates Matrix */
155 Matrix operator+(const Matrix&, const Matrix&);
156 /*! \relates Matrix */
157 Matrix operator+(const Matrix&, Matrix&&);
158 /*! \relates Matrix */
159 Matrix operator+(Matrix&&, const Matrix&);
160 /*! \relates Matrix */
161 Matrix operator+(Matrix&&, Matrix&&);
162 /*! \relates Matrix */
163 Matrix operator-(const Matrix&);
164 /*! \relates Matrix */
165 Matrix operator-(Matrix&&);
166 /*! \relates Matrix */
167 Matrix operator-(const Matrix&, const Matrix&);
168 /*! \relates Matrix */
169 Matrix operator-(const Matrix&, Matrix&&);
170 /*! \relates Matrix */
171 Matrix operator-(Matrix&&, const Matrix&);
172 /*! \relates Matrix */
173 Matrix operator-(Matrix&&, Matrix&&);
174 /*! \relates Matrix */
175 Matrix operator*(const Matrix&, Real);
176 /*! \relates Matrix */
177 Matrix operator*(Matrix&&, Real);
178 /*! \relates Matrix */
179 Matrix operator*(Real, const Matrix&);
180 /*! \relates Matrix */
181 Matrix operator*(Real, Matrix&&);
182 /*! \relates Matrix */
183 Matrix operator/(const Matrix&, Real);
184 /*! \relates Matrix */
185 Matrix operator/(Matrix&&, Real);
186
187 // vectorial products
188
189 /*! \relates Matrix */
190 Array operator*(const Array&, const Matrix&);
191 /*! \relates Matrix */
192 Array operator*(const Matrix&, const Array&);
193 /*! \relates Matrix */
194 Matrix operator*(const Matrix&, const Matrix&);
195
196 // misc. operations
197
198 /*! \relates Matrix */
199 Matrix transpose(const Matrix&);
200
201 /*! \relates Matrix */
202 Matrix outerProduct(const Array& v1, const Array& v2);
203
204 /*! \relates Matrix */
205 template <class Iterator1, class Iterator2>
206 Matrix outerProduct(Iterator1 v1begin, Iterator1 v1end, Iterator2 v2begin, Iterator2 v2end);
207
208 /*! \relates Matrix */
209 void swap(Matrix&, Matrix&) noexcept;
210
211 /*! \relates Matrix */
212 std::ostream& operator<<(std::ostream&, const Matrix&);
213
214 /*! \relates Matrix */
215 Matrix inverse(const Matrix& m);
216
217 /*! \relates Matrix */
218 Real determinant(const Matrix& m);
219
220 // inline definitions
221
222 inline Matrix::Matrix() : data_((Real*)nullptr) {}
223
224 inline Matrix::Matrix(Size rows, Size columns)
225 : data_(rows * columns > 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows),
226 columns_(columns) {}
227
228 inline Matrix::Matrix(Size rows, Size columns, Real value)
229 : data_(rows * columns > 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows),
230 columns_(columns) {
231 std::fill(first: begin(),last: end(),value: value);
232 }
233
234 template <class Iterator>
235 inline Matrix::Matrix(Size rows, Size columns, Iterator begin, Iterator end)
236 : data_(rows * columns > 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows),
237 columns_(columns) {
238 std::copy(begin, end, this->begin());
239 }
240
241 inline Matrix::Matrix(const Matrix& from)
242 : data_(!from.empty() ? new Real[from.rows_ * from.columns_] : (Real*)nullptr),
243 rows_(from.rows_), columns_(from.columns_) {
244 #if defined(QL_PATCH_MSVC) && defined(QL_DEBUG)
245 if (!from.empty())
246 #endif
247 std::copy(first: from.begin(),last: from.end(),result: begin());
248 }
249
250 inline Matrix::Matrix(Matrix&& from) noexcept
251 : data_((Real*)nullptr) {
252 swap(from);
253 }
254
255 inline Matrix::Matrix(std::initializer_list<std::initializer_list<Real>> data)
256 : data_(data.size() == 0 || data.begin()->size() == 0 ?
257 (Real*)nullptr : new Real[data.size() * data.begin()->size()]),
258 rows_(data.size()), columns_(data.size() == 0 ? 0 : data.begin()->size()) {
259 Size i=0;
260 for (const auto& row : data) {
261 #if defined(QL_EXTRA_SAFETY_CHECKS)
262 QL_REQUIRE(row.size() == columns_,
263 "a matrix needs the same number of elements for each row");
264 #endif
265 std::copy(first: row.begin(), last: row.end(), result: row_begin(i));
266 ++i;
267 }
268 }
269
270 inline Matrix& Matrix::operator=(const Matrix& from) {
271 // strong guarantee
272 Matrix temp(from);
273 swap(temp);
274 return *this;
275 }
276
277 inline Matrix& Matrix::operator=(Matrix&& from) noexcept {
278 swap(from);
279 return *this;
280 }
281
282 inline bool Matrix::operator==(const Matrix& to) const {
283 return rows_ == to.rows_ && columns_ == to.columns_ &&
284 std::equal(first1: begin(), last1: end(), first2: to.begin());
285 }
286
287 inline bool Matrix::operator!=(const Matrix& to) const {
288 return !this->operator==(to);
289 }
290
291 inline void Matrix::swap(Matrix& from) noexcept {
292 data_.swap(u&: from.data_);
293 std::swap(a&: rows_, b&: from.rows_);
294 std::swap(a&: columns_, b&: from.columns_);
295 }
296
297 inline const Matrix& Matrix::operator+=(const Matrix& m) {
298 QL_REQUIRE(rows_ == m.rows_ && columns_ == m.columns_,
299 "matrices with different sizes (" <<
300 m.rows_ << "x" << m.columns_ << ", " <<
301 rows_ << "x" << columns_ << ") cannot be "
302 "added");
303 std::transform(first1: begin(), last1: end(), first2: m.begin(), result: begin(), binary_op: std::plus<>());
304 return *this;
305 }
306
307 inline const Matrix& Matrix::operator-=(const Matrix& m) {
308 QL_REQUIRE(rows_ == m.rows_ && columns_ == m.columns_,
309 "matrices with different sizes (" <<
310 m.rows_ << "x" << m.columns_ << ", " <<
311 rows_ << "x" << columns_ << ") cannot be "
312 "subtracted");
313 std::transform(first1: begin(), last1: end(), first2: m.begin(), result: begin(), binary_op: std::minus<>());
314 return *this;
315 }
316
317 inline const Matrix& Matrix::operator*=(Real x) {
318 std::transform(first: begin(), last: end(), result: begin(), unary_op: [=](Real y) -> Real { return y * x; });
319 return *this;
320 }
321
322 inline const Matrix& Matrix::operator/=(Real x) {
323 std::transform(first: begin(),last: end(),result: begin(), unary_op: [=](Real y) -> Real { return y / x; });
324 return *this;
325 }
326
327 inline Matrix::const_iterator Matrix::begin() const {
328 return data_.get();
329 }
330
331 inline Matrix::iterator Matrix::begin() {
332 return data_.get();
333 }
334
335 inline Matrix::const_iterator Matrix::end() const {
336 return data_.get()+rows_*columns_;
337 }
338
339 inline Matrix::iterator Matrix::end() {
340 return data_.get()+rows_*columns_;
341 }
342
343 inline Matrix::const_reverse_iterator Matrix::rbegin() const {
344 return const_reverse_iterator(end());
345 }
346
347 inline Matrix::reverse_iterator Matrix::rbegin() {
348 return reverse_iterator(end());
349 }
350
351 inline Matrix::const_reverse_iterator Matrix::rend() const {
352 return const_reverse_iterator(begin());
353 }
354
355 inline Matrix::reverse_iterator Matrix::rend() {
356 return reverse_iterator(begin());
357 }
358
359 inline Matrix::const_row_iterator
360 Matrix::row_begin(Size i) const {
361 #if defined(QL_EXTRA_SAFETY_CHECKS)
362 QL_REQUIRE(i<rows_,
363 "row index (" << i << ") must be less than " << rows_ <<
364 ": matrix cannot be accessed out of range");
365 #endif
366 return data_.get()+columns_*i;
367 }
368
369 inline Matrix::row_iterator Matrix::row_begin(Size i) {
370 #if defined(QL_EXTRA_SAFETY_CHECKS)
371 QL_REQUIRE(i<rows_,
372 "row index (" << i << ") must be less than " << rows_ <<
373 ": matrix cannot be accessed out of range");
374 #endif
375 return data_.get()+columns_*i;
376 }
377
378 inline Matrix::const_row_iterator Matrix::row_end(Size i) const{
379 #if defined(QL_EXTRA_SAFETY_CHECKS)
380 QL_REQUIRE(i<rows_,
381 "row index (" << i << ") must be less than " << rows_ <<
382 ": matrix cannot be accessed out of range");
383 #endif
384 return data_.get()+columns_*(i+1);
385 }
386
387 inline Matrix::row_iterator Matrix::row_end(Size i) {
388 #if defined(QL_EXTRA_SAFETY_CHECKS)
389 QL_REQUIRE(i<rows_,
390 "row index (" << i << ") must be less than " << rows_ <<
391 ": matrix cannot be accessed out of range");
392 #endif
393 return data_.get()+columns_*(i+1);
394 }
395
396 inline Matrix::const_reverse_row_iterator
397 Matrix::row_rbegin(Size i) const {
398 return const_reverse_row_iterator(row_end(i));
399 }
400
401 inline Matrix::reverse_row_iterator Matrix::row_rbegin(Size i) {
402 return reverse_row_iterator(row_end(i));
403 }
404
405 inline Matrix::const_reverse_row_iterator
406 Matrix::row_rend(Size i) const {
407 return const_reverse_row_iterator(row_begin(i));
408 }
409
410 inline Matrix::reverse_row_iterator Matrix::row_rend(Size i) {
411 return reverse_row_iterator(row_begin(i));
412 }
413
414 inline Matrix::const_column_iterator
415 Matrix::column_begin(Size i) const {
416 #if defined(QL_EXTRA_SAFETY_CHECKS)
417 QL_REQUIRE(i<columns_,
418 "column index (" << i << ") must be less than " << columns_ <<
419 ": matrix cannot be accessed out of range");
420 #endif
421 return const_column_iterator(data_.get()+i,columns_);
422 }
423
424 inline Matrix::column_iterator Matrix::column_begin(Size i) {
425 #if defined(QL_EXTRA_SAFETY_CHECKS)
426 QL_REQUIRE(i<columns_,
427 "column index (" << i << ") must be less than " << columns_ <<
428 ": matrix cannot be accessed out of range");
429 #endif
430 return column_iterator(data_.get()+i,columns_);
431 }
432
433 inline Matrix::const_column_iterator
434 Matrix::column_end(Size i) const {
435 #if defined(QL_EXTRA_SAFETY_CHECKS)
436 QL_REQUIRE(i<columns_,
437 "column index (" << i << ") must be less than " << columns_ <<
438 ": matrix cannot be accessed out of range");
439 #endif
440 return const_column_iterator(data_.get()+i+rows_*columns_,columns_);
441 }
442
443 inline Matrix::column_iterator Matrix::column_end(Size i) {
444 #if defined(QL_EXTRA_SAFETY_CHECKS)
445 QL_REQUIRE(i<columns_,
446 "column index (" << i << ") must be less than " << columns_ <<
447 ": matrix cannot be accessed out of range");
448 #endif
449 return column_iterator(data_.get()+i+rows_*columns_,columns_);
450 }
451
452 inline Matrix::const_reverse_column_iterator
453 Matrix::column_rbegin(Size i) const {
454 return const_reverse_column_iterator(column_end(i));
455 }
456
457 inline Matrix::reverse_column_iterator
458 Matrix::column_rbegin(Size i) {
459 return reverse_column_iterator(column_end(i));
460 }
461
462 inline Matrix::const_reverse_column_iterator
463 Matrix::column_rend(Size i) const {
464 return const_reverse_column_iterator(column_begin(i));
465 }
466
467 inline Matrix::reverse_column_iterator
468 Matrix::column_rend(Size i) {
469 return reverse_column_iterator(column_begin(i));
470 }
471
472 inline Matrix::const_row_iterator
473 Matrix::operator[](Size i) const {
474 return row_begin(i);
475 }
476
477 inline Matrix::const_row_iterator
478 Matrix::at(Size i) const {
479 QL_REQUIRE(i < rows_, "matrix access out of range");
480 return row_begin(i);
481 }
482
483 inline Matrix::row_iterator Matrix::operator[](Size i) {
484 return row_begin(i);
485 }
486
487 inline Matrix::row_iterator Matrix::at(Size i) {
488 QL_REQUIRE(i < rows_, "matrix access out of range");
489 return row_begin(i);
490 }
491
492 inline Array Matrix::diagonal() const {
493 Size arraySize = std::min<Size>(a: rows(), b: columns());
494 Array tmp(arraySize);
495 for(Size i = 0; i < arraySize; i++)
496 tmp[i] = (*this)[i][i];
497 return tmp;
498 }
499
500 inline Real &Matrix::operator()(Size i, Size j) const {
501 return data_[i*columns()+j];
502 }
503
504 inline Size Matrix::rows() const {
505 return rows_;
506 }
507
508 inline Size Matrix::columns() const {
509 return columns_;
510 }
511
512 inline Size Matrix::size1() const {
513 return rows();
514 }
515
516 inline Size Matrix::size2() const {
517 return columns();
518 }
519
520 inline bool Matrix::empty() const {
521 return rows_ == 0 || columns_ == 0;
522 }
523
524 inline Matrix operator+(const Matrix& m1, const Matrix& m2) {
525 QL_REQUIRE(m1.rows() == m2.rows() &&
526 m1.columns() == m2.columns(),
527 "matrices with different sizes (" <<
528 m1.rows() << "x" << m1.columns() << ", " <<
529 m2.rows() << "x" << m2.columns() << ") cannot be "
530 "added");
531 Matrix temp(m1.rows(),m1.columns());
532 std::transform(first1: m1.begin(), last1: m1.end(), first2: m2.begin(), result: temp.begin(), binary_op: std::plus<>());
533 return temp;
534 }
535
536 inline Matrix operator+(const Matrix& m1, Matrix&& m2) {
537 QL_REQUIRE(m1.rows() == m2.rows() &&
538 m1.columns() == m2.columns(),
539 "matrices with different sizes (" <<
540 m1.rows() << "x" << m1.columns() << ", " <<
541 m2.rows() << "x" << m2.columns() << ") cannot be "
542 "added");
543 std::transform(first1: m1.begin(), last1: m1.end(), first2: m2.begin(), result: m2.begin(), binary_op: std::plus<>());
544 return std::move(m2);
545 }
546
547 inline Matrix operator+(Matrix&& m1, const Matrix& m2) {
548 QL_REQUIRE(m1.rows() == m2.rows() &&
549 m1.columns() == m2.columns(),
550 "matrices with different sizes (" <<
551 m1.rows() << "x" << m1.columns() << ", " <<
552 m2.rows() << "x" << m2.columns() << ") cannot be "
553 "added");
554 std::transform(first1: m1.begin(), last1: m1.end(), first2: m2.begin(), result: m1.begin(), binary_op: std::plus<>());
555 return std::move(m1);
556 }
557
558 inline Matrix operator+(Matrix&& m1, Matrix&& m2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved)
559 QL_REQUIRE(m1.rows() == m2.rows() &&
560 m1.columns() == m2.columns(),
561 "matrices with different sizes (" <<
562 m1.rows() << "x" << m1.columns() << ", " <<
563 m2.rows() << "x" << m2.columns() << ") cannot be "
564 "added");
565 std::transform(first1: m1.begin(), last1: m1.end(), first2: m2.begin(), result: m1.begin(), binary_op: std::plus<>());
566 return std::move(m1);
567 }
568
569 inline Matrix operator-(const Matrix& m1) {
570 Matrix temp(m1.rows(), m1.columns());
571 std::transform(first: m1.begin(), last: m1.end(), result: temp.begin(), unary_op: std::negate<>());
572 return temp;
573 }
574
575 inline Matrix operator-(Matrix&& m1) {
576 std::transform(first: m1.begin(), last: m1.end(), result: m1.begin(), unary_op: std::negate<>());
577 return std::move(m1);
578 }
579
580 inline Matrix operator-(const Matrix& m1, const Matrix& m2) {
581 QL_REQUIRE(m1.rows() == m2.rows() &&
582 m1.columns() == m2.columns(),
583 "matrices with different sizes (" <<
584 m1.rows() << "x" << m1.columns() << ", " <<
585 m2.rows() << "x" << m2.columns() << ") cannot be "
586 "subtracted");
587 Matrix temp(m1.rows(),m1.columns());
588 std::transform(first1: m1.begin(), last1: m1.end(), first2: m2.begin(), result: temp.begin(), binary_op: std::minus<>());
589 return temp;
590 }
591
592 inline Matrix operator-(const Matrix& m1, Matrix&& m2) {
593 QL_REQUIRE(m1.rows() == m2.rows() &&
594 m1.columns() == m2.columns(),
595 "matrices with different sizes (" <<
596 m1.rows() << "x" << m1.columns() << ", " <<
597 m2.rows() << "x" << m2.columns() << ") cannot be "
598 "subtracted");
599 std::transform(first1: m1.begin(), last1: m1.end(), first2: m2.begin(), result: m2.begin(), binary_op: std::minus<>());
600 return std::move(m2);
601 }
602
603 inline Matrix operator-(Matrix&& m1, const Matrix& m2) {
604 QL_REQUIRE(m1.rows() == m2.rows() &&
605 m1.columns() == m2.columns(),
606 "matrices with different sizes (" <<
607 m1.rows() << "x" << m1.columns() << ", " <<
608 m2.rows() << "x" << m2.columns() << ") cannot be "
609 "subtracted");
610 std::transform(first1: m1.begin(), last1: m1.end(), first2: m2.begin(), result: m1.begin(), binary_op: std::minus<>());
611 return std::move(m1);
612 }
613
614 inline Matrix operator-(Matrix&& m1, Matrix&& m2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved)
615 QL_REQUIRE(m1.rows() == m2.rows() &&
616 m1.columns() == m2.columns(),
617 "matrices with different sizes (" <<
618 m1.rows() << "x" << m1.columns() << ", " <<
619 m2.rows() << "x" << m2.columns() << ") cannot be "
620 "subtracted");
621 std::transform(first1: m1.begin(), last1: m1.end(), first2: m2.begin(), result: m1.begin(), binary_op: std::minus<>());
622 return std::move(m1);
623 }
624
625 inline Matrix operator*(const Matrix& m, Real x) {
626 Matrix temp(m.rows(),m.columns());
627 std::transform(first: m.begin(), last: m.end(), result: temp.begin(), unary_op: [=](Real y) -> Real { return y * x; });
628 return temp;
629 }
630
631 inline Matrix operator*(Matrix&& m, Real x) {
632 std::transform(first: m.begin(), last: m.end(), result: m.begin(), unary_op: [=](Real y) -> Real { return y * x; });
633 return std::move(m);
634 }
635
636 inline Matrix operator*(Real x, const Matrix& m) {
637 Matrix temp(m.rows(),m.columns());
638 std::transform(first: m.begin(), last: m.end(), result: temp.begin(), unary_op: [=](Real y) -> Real { return x * y; });
639 return temp;
640 }
641
642 inline Matrix operator*(Real x, Matrix&& m) {
643 std::transform(first: m.begin(), last: m.end(), result: m.begin(), unary_op: [=](Real y) -> Real { return x * y; });
644 return std::move(m);
645 }
646
647 inline Matrix operator/(const Matrix& m, Real x) {
648 Matrix temp(m.rows(),m.columns());
649 std::transform(first: m.begin(), last: m.end(), result: temp.begin(), unary_op: [=](Real y) -> Real { return y / x; });
650 return temp;
651 }
652
653 inline Matrix operator/(Matrix&& m, Real x) {
654 std::transform(first: m.begin(), last: m.end(), result: m.begin(), unary_op: [=](Real y) -> Real { return y / x; });
655 return std::move(m);
656 }
657
658 inline Array operator*(const Array& v, const Matrix& m) {
659 QL_REQUIRE(v.size() == m.rows(),
660 "vectors and matrices with different sizes ("
661 << v.size() << ", " << m.rows() << "x" << m.columns() <<
662 ") cannot be multiplied");
663 Array result(m.columns());
664 for (Size i=0; i<result.size(); i++)
665 result[i] =
666 std::inner_product(first1: v.begin(),last1: v.end(),
667 first2: m.column_begin(i),init: Real(0.0));
668 return result;
669 }
670
671 inline Array operator*(const Matrix& m, const Array& v) {
672 QL_REQUIRE(v.size() == m.columns(),
673 "vectors and matrices with different sizes ("
674 << v.size() << ", " << m.rows() << "x" << m.columns() <<
675 ") cannot be multiplied");
676 Array result(m.rows());
677 for (Size i=0; i<result.size(); i++)
678 result[i] =
679 std::inner_product(first1: v.begin(),last1: v.end(),first2: m.row_begin(i),init: Real(0.0));
680 return result;
681 }
682
683 inline Matrix operator*(const Matrix& m1, const Matrix& m2) {
684 QL_REQUIRE(m1.columns() == m2.rows(),
685 "matrices with different sizes (" <<
686 m1.rows() << "x" << m1.columns() << ", " <<
687 m2.rows() << "x" << m2.columns() << ") cannot be "
688 "multiplied");
689 Matrix result(m1.rows(),m2.columns(),0.0);
690 for (Size i=0; i<result.rows(); ++i) {
691 for (Size k=0; k<m1.columns(); ++k) {
692 for (Size j=0; j<result.columns(); ++j) {
693 result[i][j] += m1[i][k]*m2[k][j];
694 }
695 }
696 }
697 return result;
698 }
699
700 inline Matrix transpose(const Matrix& m) {
701 Matrix result(m.columns(),m.rows());
702 #if defined(QL_PATCH_MSVC) && defined(QL_DEBUG)
703 if (!m.empty())
704 #endif
705 for (Size i=0; i<m.rows(); i++)
706 std::copy(first: m.row_begin(i),last: m.row_end(i),result: result.column_begin(i));
707 return result;
708 }
709
710 inline Matrix outerProduct(const Array& v1, const Array& v2) {
711 return outerProduct(v1begin: v1.begin(), v1end: v1.end(), v2begin: v2.begin(), v2end: v2.end());
712 }
713
714 template <class Iterator1, class Iterator2>
715 inline Matrix outerProduct(Iterator1 v1begin, Iterator1 v1end, Iterator2 v2begin, Iterator2 v2end) {
716
717 Size size1 = std::distance(v1begin, v1end);
718 QL_REQUIRE(size1>0, "null first vector");
719
720 Size size2 = std::distance(v2begin, v2end);
721 QL_REQUIRE(size2>0, "null second vector");
722
723 Matrix result(size1, size2);
724
725 for (Size i=0; v1begin!=v1end; i++, v1begin++)
726 std::transform(v2begin, v2end, result.row_begin(i),
727 [=](Real y) -> Real { return y * (*v1begin); });
728
729 return result;
730 }
731
732 inline void swap(Matrix& m1, Matrix& m2) noexcept {
733 m1.swap(from&: m2);
734 }
735
736 inline std::ostream& operator<<(std::ostream& out, const Matrix& m) {
737 std::streamsize width = out.width();
738 for (Size i=0; i<m.rows(); i++) {
739 out << "| ";
740 for (Size j=0; j<m.columns(); j++)
741 out << std::setw(int(width)) << m[i][j] << " ";
742 out << "|\n";
743 }
744 return out;
745 }
746
747}
748
749
750#endif
751

source code of quantlib/ql/math/matrix.hpp