42 |
26 Feb 04 |
jari |
// $Id$ |
12 |
19 Jun 03 |
daniel |
2 |
|
570 |
05 Apr 06 |
jari |
3 |
/* |
570 |
05 Apr 06 |
jari |
Copyright (C) 2003 Daniel Dalevi, Peter Johansson |
2119 |
12 Dec 09 |
peter |
Copyright (C) 2004 Jari Häkkinen, Peter Johansson |
2119 |
12 Dec 09 |
peter |
Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér |
2119 |
12 Dec 09 |
peter |
Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson |
4359 |
23 Aug 23 |
peter |
Copyright (C) 2009, 2010, 2011, 2012, 2017, 2021, 2022, 2023 Peter Johansson |
570 |
05 Apr 06 |
jari |
9 |
|
1437 |
25 Aug 08 |
peter |
This file is part of the yat library, http://dev.thep.lu.se/yat |
570 |
05 Apr 06 |
jari |
11 |
|
675 |
10 Oct 06 |
jari |
The yat library is free software; you can redistribute it and/or |
675 |
10 Oct 06 |
jari |
modify it under the terms of the GNU General Public License as |
1486 |
09 Sep 08 |
jari |
published by the Free Software Foundation; either version 3 of the |
675 |
10 Oct 06 |
jari |
License, or (at your option) any later version. |
570 |
05 Apr 06 |
jari |
16 |
|
675 |
10 Oct 06 |
jari |
The yat library is distributed in the hope that it will be useful, |
675 |
10 Oct 06 |
jari |
but WITHOUT ANY WARRANTY; without even the implied warranty of |
675 |
10 Oct 06 |
jari |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
570 |
05 Apr 06 |
jari |
General Public License for more details. |
570 |
05 Apr 06 |
jari |
21 |
|
570 |
05 Apr 06 |
jari |
You should have received a copy of the GNU General Public License |
1487 |
10 Sep 08 |
jari |
along with yat. If not, see <http://www.gnu.org/licenses/>. |
570 |
05 Apr 06 |
jari |
24 |
*/ |
570 |
05 Apr 06 |
jari |
25 |
|
2881 |
18 Nov 12 |
peter |
26 |
#include <config.h> |
2881 |
18 Nov 12 |
peter |
27 |
|
1121 |
22 Feb 08 |
peter |
28 |
#include "Matrix.h" |
2210 |
05 Mar 10 |
peter |
29 |
|
3657 |
13 Jul 17 |
peter |
30 |
#include "DiagonalMatrix.h" |
2210 |
05 Mar 10 |
peter |
31 |
#include "Exception.h" |
3657 |
13 Jul 17 |
peter |
32 |
#include "SVD.h" |
1121 |
22 Feb 08 |
peter |
33 |
#include "Vector.h" |
1017 |
01 Feb 08 |
peter |
34 |
#include "VectorBase.h" |
1027 |
02 Feb 08 |
peter |
35 |
#include "VectorConstView.h" |
1015 |
01 Feb 08 |
peter |
36 |
#include "VectorView.h" |
680 |
11 Oct 06 |
jari |
37 |
#include "utility.h" |
273 |
14 Apr 05 |
peter |
38 |
|
1792 |
11 Feb 09 |
peter |
39 |
#include <gsl/gsl_blas.h> |
1792 |
11 Feb 09 |
peter |
40 |
|
2325 |
22 Sep 10 |
peter |
41 |
#include <algorithm> |
792 |
11 Mar 07 |
jari |
42 |
#include <cassert> |
1783 |
08 Feb 09 |
peter |
43 |
#include <climits> |
293 |
26 Apr 05 |
peter |
44 |
#include <cmath> |
1792 |
11 Feb 09 |
peter |
45 |
#include <iostream> |
1920 |
24 Apr 09 |
peter |
46 |
#include <limits> |
42 |
26 Feb 04 |
jari |
47 |
#include <sstream> |
17 |
08 Jul 03 |
peter |
48 |
#include <vector> |
17 |
08 Jul 03 |
peter |
49 |
|
42 |
26 Feb 04 |
jari |
50 |
namespace theplu { |
680 |
11 Oct 06 |
jari |
51 |
namespace yat { |
616 |
31 Aug 06 |
jari |
52 |
namespace utility { |
12 |
19 Jun 03 |
daniel |
53 |
|
12 |
19 Jun 03 |
daniel |
54 |
|
1121 |
22 Feb 08 |
peter |
55 |
Matrix::Matrix(void) |
4129 |
19 Jan 22 |
peter |
56 |
: m_(nullptr) |
703 |
18 Dec 06 |
jari |
57 |
{ |
703 |
18 Dec 06 |
jari |
58 |
} |
12 |
19 Jun 03 |
daniel |
59 |
|
703 |
18 Dec 06 |
jari |
60 |
|
1121 |
22 Feb 08 |
peter |
61 |
Matrix::Matrix(const size_t& r, const size_t& c, double init_value) |
4129 |
19 Jan 22 |
peter |
62 |
: m_(nullptr) |
703 |
18 Dec 06 |
jari |
63 |
{ |
4124 |
13 Jan 22 |
peter |
64 |
assert(static_cast<bool>(r) == static_cast<bool>(c)); |
1598 |
24 Oct 08 |
peter |
65 |
resize(r,c,init_value); |
4124 |
13 Jan 22 |
peter |
66 |
assert(rows() == r); |
4124 |
13 Jan 22 |
peter |
67 |
assert(columns() == c); |
703 |
18 Dec 06 |
jari |
68 |
} |
703 |
18 Dec 06 |
jari |
69 |
|
703 |
18 Dec 06 |
jari |
70 |
|
1121 |
22 Feb 08 |
peter |
71 |
Matrix::Matrix(const Matrix& o) |
4129 |
19 Jan 22 |
peter |
72 |
: m_(nullptr) |
703 |
18 Dec 06 |
jari |
73 |
{ |
4124 |
13 Jan 22 |
peter |
74 |
m_ = detail::create_copy(o.gsl_matrix_p()); |
703 |
18 Dec 06 |
jari |
75 |
} |
703 |
18 Dec 06 |
jari |
76 |
|
703 |
18 Dec 06 |
jari |
77 |
|
4124 |
13 Jan 22 |
peter |
78 |
Matrix::Matrix(const MatrixBase& o) |
4129 |
19 Jan 22 |
peter |
79 |
: m_(nullptr) |
4124 |
13 Jan 22 |
peter |
80 |
{ |
4124 |
13 Jan 22 |
peter |
81 |
m_ = detail::create_copy(o.gsl_matrix_p()); |
4124 |
13 Jan 22 |
peter |
82 |
} |
4124 |
13 Jan 22 |
peter |
83 |
|
4124 |
13 Jan 22 |
peter |
84 |
|
3691 |
14 Sep 17 |
peter |
85 |
Matrix::Matrix(Matrix&& other) noexcept |
4129 |
19 Jan 22 |
peter |
86 |
: m_(nullptr) |
3585 |
19 Jan 17 |
peter |
87 |
{ |
3585 |
19 Jan 17 |
peter |
88 |
std::swap(m_, other.m_); |
3585 |
19 Jan 17 |
peter |
89 |
} |
3585 |
19 Jan 17 |
peter |
90 |
|
3585 |
19 Jan 17 |
peter |
91 |
|
4124 |
13 Jan 22 |
peter |
92 |
Matrix::Matrix(MatrixMutable&& other) |
4129 |
19 Jan 22 |
peter |
93 |
: m_(nullptr) |
4124 |
13 Jan 22 |
peter |
94 |
{ |
4378 |
11 Oct 23 |
peter |
95 |
move_assign(std::move(other)); |
4124 |
13 Jan 22 |
peter |
96 |
} |
4124 |
13 Jan 22 |
peter |
97 |
|
4124 |
13 Jan 22 |
peter |
98 |
|
42 |
26 Feb 04 |
jari |
// Constructor that gets data from istream |
2689 |
27 Feb 12 |
peter |
100 |
Matrix::Matrix(std::istream& is, char sep) |
4129 |
19 Jan 22 |
peter |
101 |
: m_(nullptr) |
42 |
26 Feb 04 |
jari |
102 |
{ |
1147 |
25 Feb 08 |
peter |
103 |
if (!is.good()) |
1147 |
25 Feb 08 |
peter |
104 |
throw utility::IO_error("Matrix: istream is not good"); |
1147 |
25 Feb 08 |
peter |
105 |
|
42 |
26 Feb 04 |
jari |
// read the data file and store in stl vectors (dynamically |
42 |
26 Feb 04 |
jari |
// expandable) |
42 |
26 Feb 04 |
jari |
108 |
std::vector<std::vector<double> > data_matrix; |
1927 |
30 Apr 09 |
peter |
109 |
try { |
1928 |
30 Apr 09 |
peter |
110 |
load(is, data_matrix, sep, '\n', true); |
1927 |
30 Apr 09 |
peter |
111 |
} |
1928 |
30 Apr 09 |
peter |
112 |
catch (utility::IO_error& e) { |
4179 |
02 Jun 22 |
peter |
113 |
std::stringstream ss; |
4179 |
02 Jun 22 |
peter |
114 |
ss << "Matrix(std::istream&): invalid dimensions"; |
4060 |
10 May 21 |
peter |
115 |
std::throw_with_nested(IO_error(ss.str())); |
1928 |
30 Apr 09 |
peter |
116 |
} |
2210 |
05 Mar 10 |
peter |
117 |
catch (runtime_error& e) { |
4179 |
02 Jun 22 |
peter |
118 |
std::stringstream ss; |
4179 |
02 Jun 22 |
peter |
119 |
ss << "Matrix(std::istream&): invalid matrix element"; |
4060 |
10 May 21 |
peter |
120 |
std::throw_with_nested(IO_error(ss.str())); |
1927 |
30 Apr 09 |
peter |
121 |
} |
1927 |
30 Apr 09 |
peter |
122 |
|
1928 |
30 Apr 09 |
peter |
123 |
unsigned int nof_rows = data_matrix.size(); |
1147 |
25 Feb 08 |
peter |
// if stream was empty, create nothing |
1928 |
30 Apr 09 |
peter |
125 |
if (!nof_rows) |
1147 |
25 Feb 08 |
peter |
126 |
return; |
1147 |
25 Feb 08 |
peter |
127 |
|
1928 |
30 Apr 09 |
peter |
128 |
unsigned int nof_columns=data_matrix[0].size(); |
1928 |
30 Apr 09 |
peter |
129 |
|
42 |
26 Feb 04 |
jari |
// convert the data to a gsl matrix |
3651 |
28 Jun 17 |
peter |
131 |
m_ = detail::allocate_matrix(nof_rows, nof_columns); |
4124 |
13 Jan 22 |
peter |
132 |
assert(columns() == nof_columns); |
1927 |
30 Apr 09 |
peter |
133 |
for (size_t i=0; i<data_matrix.size(); ++i) { |
1927 |
30 Apr 09 |
peter |
134 |
assert(data_matrix[i].size()==columns()); |
1928 |
30 Apr 09 |
peter |
135 |
assert(i<rows()); |
1928 |
30 Apr 09 |
peter |
136 |
std::copy(data_matrix[i].begin(), data_matrix[i].end(), begin_row(i)); |
1927 |
30 Apr 09 |
peter |
137 |
} |
42 |
26 Feb 04 |
jari |
138 |
} |
12 |
19 Jun 03 |
daniel |
139 |
|
12 |
19 Jun 03 |
daniel |
140 |
|
4174 |
18 May 22 |
peter |
141 |
Matrix::~Matrix(void) |
4174 |
18 May 22 |
peter |
142 |
{ |
4174 |
18 May 22 |
peter |
143 |
detail::deallocate(m_); |
4174 |
18 May 22 |
peter |
144 |
} |
4174 |
18 May 22 |
peter |
145 |
|
4174 |
18 May 22 |
peter |
146 |
|
4129 |
19 Jan 22 |
peter |
147 |
void Matrix::copy_assign(const gsl_matrix* rhs) |
4129 |
19 Jan 22 |
peter |
148 |
{ |
4129 |
19 Jan 22 |
peter |
// copy via blas_result and then swap result to m_ |
4129 |
19 Jan 22 |
peter |
150 |
detail::copy(blas_result_, rhs); |
4129 |
19 Jan 22 |
peter |
151 |
std::swap(blas_result_, m_); |
4129 |
19 Jan 22 |
peter |
152 |
} |
4129 |
19 Jan 22 |
peter |
153 |
|
4129 |
19 Jan 22 |
peter |
154 |
|
4129 |
19 Jan 22 |
peter |
155 |
gsl_matrix* Matrix::gsl_matrix_p(void) |
4129 |
19 Jan 22 |
peter |
156 |
{ |
4129 |
19 Jan 22 |
peter |
157 |
return m_; |
4129 |
19 Jan 22 |
peter |
158 |
} |
4129 |
19 Jan 22 |
peter |
159 |
|
4129 |
19 Jan 22 |
peter |
160 |
|
4129 |
19 Jan 22 |
peter |
161 |
const gsl_matrix* Matrix::gsl_matrix_p(void) const |
4129 |
19 Jan 22 |
peter |
162 |
{ |
4129 |
19 Jan 22 |
peter |
163 |
return m_; |
4129 |
19 Jan 22 |
peter |
164 |
} |
4129 |
19 Jan 22 |
peter |
165 |
|
4129 |
19 Jan 22 |
peter |
166 |
|
4129 |
19 Jan 22 |
peter |
167 |
void Matrix::move_assign(MatrixMutable&& rhs) |
4129 |
19 Jan 22 |
peter |
168 |
{ |
4129 |
19 Jan 22 |
peter |
169 |
detail::Mover mover(*this); |
4129 |
19 Jan 22 |
peter |
170 |
mover(rhs); |
4129 |
19 Jan 22 |
peter |
171 |
} |
4129 |
19 Jan 22 |
peter |
172 |
|
4129 |
19 Jan 22 |
peter |
173 |
|
4129 |
19 Jan 22 |
peter |
174 |
void Matrix::move_assign(gsl_matrix*&& rhs) |
4129 |
19 Jan 22 |
peter |
175 |
{ |
4129 |
19 Jan 22 |
peter |
// if there is a use case for when rhs is a view, we could allow |
4129 |
19 Jan 22 |
peter |
// that with an if statement, but ATM it seems a waste. |
4129 |
19 Jan 22 |
peter |
178 |
assert(!rhs || rhs->owner); |
4129 |
19 Jan 22 |
peter |
179 |
std::swap(m_, rhs); |
4129 |
19 Jan 22 |
peter |
180 |
} |
4129 |
19 Jan 22 |
peter |
181 |
|
4129 |
19 Jan 22 |
peter |
182 |
|
1121 |
22 Feb 08 |
peter |
183 |
void Matrix::resize(size_t r, size_t c, double init_value) |
808 |
15 Mar 07 |
peter |
184 |
{ |
3651 |
28 Jun 17 |
peter |
185 |
detail::reallocate(m_, r, c); |
808 |
15 Mar 07 |
peter |
186 |
|
3651 |
28 Jun 17 |
peter |
187 |
if (m_) |
1172 |
27 Feb 08 |
peter |
188 |
all(init_value); |
808 |
15 Mar 07 |
peter |
189 |
|
810 |
15 Mar 07 |
jari |
// no need to delete blas_result_ if the number of rows fit, it |
810 |
15 Mar 07 |
jari |
// may be useful later. |
4124 |
13 Jan 22 |
peter |
192 |
|
4124 |
13 Jan 22 |
peter |
// Peter, actually it might be useful regardless of dimension |
4124 |
13 Jan 22 |
peter |
// since if we assign from MatrixExpression the rhs might have any |
4124 |
13 Jan 22 |
peter |
// dimension. The condition above was based on the fact that |
4124 |
13 Jan 22 |
peter |
// blas_result_ was mostly used in operator*= where number of rows |
4124 |
13 Jan 22 |
peter |
// is constant. |
810 |
15 Mar 07 |
jari |
198 |
if (blas_result_ && (blas_result_->size1!=rows())) { |
3651 |
28 Jun 17 |
peter |
199 |
detail::deallocate(blas_result_); |
810 |
15 Mar 07 |
jari |
200 |
} |
810 |
15 Mar 07 |
jari |
201 |
} |
808 |
15 Mar 07 |
peter |
202 |
|
810 |
15 Mar 07 |
jari |
203 |
|
1121 |
22 Feb 08 |
peter |
204 |
void Matrix::transpose(void) |
42 |
26 Feb 04 |
jari |
205 |
{ |
792 |
11 Mar 07 |
jari |
206 |
assert(m_); |
4124 |
13 Jan 22 |
peter |
207 |
if (!m_) // allow empty matrix |
4124 |
13 Jan 22 |
peter |
208 |
return; |
420 |
02 Dec 05 |
jari |
209 |
if (columns()==rows()) |
754 |
17 Feb 07 |
jari |
210 |
gsl_matrix_transpose(m_); // this never fails |
420 |
02 Dec 05 |
jari |
211 |
else { |
3651 |
28 Jun 17 |
peter |
212 |
gsl_matrix* transposed = detail::allocate_matrix(columns(), rows()); |
754 |
17 Feb 07 |
jari |
// next line never fails if allocation above succeeded. |
420 |
02 Dec 05 |
jari |
214 |
gsl_matrix_transpose_memcpy(transposed,m_); |
4124 |
13 Jan 22 |
peter |
215 |
detail::deallocate(m_); |
1098 |
18 Feb 08 |
jari |
216 |
m_ = transposed; |
3651 |
28 Jun 17 |
peter |
217 |
detail::deallocate(blas_result_); |
420 |
02 Dec 05 |
jari |
218 |
} |
42 |
26 Feb 04 |
jari |
219 |
} |
21 |
11 Jul 03 |
daniel |
220 |
|
12 |
19 Jun 03 |
daniel |
221 |
|
4129 |
19 Jan 22 |
peter |
222 |
void Matrix::visit(detail::Mover& mover) |
12 |
19 Jun 03 |
daniel |
223 |
{ |
4129 |
19 Jan 22 |
peter |
224 |
std::swap(mover.lhs().m_, m_); |
42 |
26 Feb 04 |
jari |
225 |
} |
12 |
19 Jun 03 |
daniel |
226 |
|
12 |
19 Jun 03 |
daniel |
227 |
|
4129 |
19 Jan 22 |
peter |
228 |
const Matrix& Matrix::operator=( const Matrix& other ) |
4124 |
13 Jan 22 |
peter |
229 |
{ |
4124 |
13 Jan 22 |
peter |
230 |
copy_assign(other); |
4124 |
13 Jan 22 |
peter |
231 |
return *this; |
4124 |
13 Jan 22 |
peter |
232 |
} |
4124 |
13 Jan 22 |
peter |
233 |
|
4124 |
13 Jan 22 |
peter |
234 |
|
3585 |
19 Jan 17 |
peter |
235 |
Matrix& Matrix::operator=(Matrix&& other) |
3585 |
19 Jan 17 |
peter |
236 |
{ |
4129 |
19 Jan 22 |
peter |
237 |
std::swap(m_, other.m_); |
3585 |
19 Jan 17 |
peter |
238 |
return *this; |
3585 |
19 Jan 17 |
peter |
239 |
} |
3585 |
19 Jan 17 |
peter |
240 |
|
3585 |
19 Jan 17 |
peter |
241 |
|
1121 |
22 Feb 08 |
peter |
242 |
const Matrix& Matrix::operator*=(const Matrix& other) |
42 |
26 Feb 04 |
jari |
243 |
{ |
792 |
11 Mar 07 |
jari |
244 |
assert(m_); |
3651 |
28 Jun 17 |
peter |
245 |
assert(other.columns()); |
4124 |
13 Jan 22 |
peter |
246 |
assert(other.rows() == columns()); |
3651 |
28 Jun 17 |
peter |
247 |
detail::reallocate(blas_result_, rows(), other.columns()); |
4124 |
13 Jan 22 |
peter |
248 |
assert(detail::rows(blas_result_) == detail::rows(m_)); |
4124 |
13 Jan 22 |
peter |
249 |
assert(detail::columns(blas_result_) == detail::columns(other.m_)); |
2661 |
21 Nov 11 |
peter |
250 |
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, |
2325 |
22 Sep 10 |
peter |
251 |
blas_result_); |
2325 |
22 Sep 10 |
peter |
252 |
std::swap(m_, blas_result_); |
420 |
02 Dec 05 |
jari |
253 |
return *this; |
42 |
26 Feb 04 |
jari |
254 |
} |
20 |
11 Jul 03 |
daniel |
255 |
|
12 |
19 Jun 03 |
daniel |
256 |
|
4140 |
29 Jan 22 |
peter |
257 |
void inverse_svd(const MatrixBase& A, Matrix& result) |
3657 |
13 Jul 17 |
peter |
258 |
{ |
3657 |
13 Jul 17 |
peter |
259 |
assert(A.rows() == A.columns()); |
3657 |
13 Jul 17 |
peter |
260 |
SVD svd(A); |
3657 |
13 Jul 17 |
peter |
261 |
svd.decompose(); |
4286 |
30 Jan 23 |
peter |
262 |
svd.inverse(result); |
3657 |
13 Jul 17 |
peter |
263 |
assert(result.columns() == result.rows()); |
3657 |
13 Jul 17 |
peter |
264 |
assert(result.rows() == A.rows()); |
3657 |
13 Jul 17 |
peter |
265 |
} |
3657 |
13 Jul 17 |
peter |
266 |
|
3657 |
13 Jul 17 |
peter |
267 |
|
4138 |
28 Jan 22 |
peter |
268 |
bool nan(const MatrixBase& templat, Matrix& flag) |
774 |
01 Mar 07 |
jari |
269 |
{ |
1786 |
09 Feb 09 |
peter |
270 |
if (templat.rows()!=flag.rows() || templat.columns()!=flag.columns()) |
1786 |
09 Feb 09 |
peter |
271 |
flag.resize(templat.rows(),templat.columns(),1.0); |
1807 |
18 Feb 09 |
peter |
272 |
return binary_weight(templat.begin(), templat.end(), flag.begin()); |
774 |
01 Mar 07 |
jari |
273 |
} |
774 |
01 Mar 07 |
jari |
274 |
|
774 |
01 Mar 07 |
jari |
275 |
|
4129 |
19 Jan 22 |
peter |
276 |
namespace detail { |
4129 |
19 Jan 22 |
peter |
277 |
|
4129 |
19 Jan 22 |
peter |
278 |
Mover::Mover(Matrix& m) |
4129 |
19 Jan 22 |
peter |
279 |
: lhs_(m) {} |
4129 |
19 Jan 22 |
peter |
280 |
|
4129 |
19 Jan 22 |
peter |
281 |
|
4129 |
19 Jan 22 |
peter |
282 |
void Mover::operator()(MatrixMutable& rhs) |
4129 |
19 Jan 22 |
peter |
283 |
{ |
4129 |
19 Jan 22 |
peter |
284 |
rhs.visit(*this); |
4129 |
19 Jan 22 |
peter |
285 |
} |
4129 |
19 Jan 22 |
peter |
286 |
|
4129 |
19 Jan 22 |
peter |
287 |
|
4129 |
19 Jan 22 |
peter |
288 |
void Mover::copy_assign(const MatrixMutable& rhs) |
4129 |
19 Jan 22 |
peter |
289 |
{ |
4129 |
19 Jan 22 |
peter |
290 |
lhs_.copy_assign(rhs); |
4129 |
19 Jan 22 |
peter |
291 |
} |
4129 |
19 Jan 22 |
peter |
292 |
|
4129 |
19 Jan 22 |
peter |
293 |
|
4129 |
19 Jan 22 |
peter |
294 |
Matrix& Mover::lhs(void) |
4129 |
19 Jan 22 |
peter |
295 |
{ |
4129 |
19 Jan 22 |
peter |
296 |
return lhs_; |
4129 |
19 Jan 22 |
peter |
297 |
} |
4129 |
19 Jan 22 |
peter |
298 |
} |
4129 |
19 Jan 22 |
peter |
299 |
|
680 |
11 Oct 06 |
jari |
300 |
}}} // of namespace utility, yat and thep |