PLaSK library

Loading...
Searching...
No Matches
matrices.hpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
17#ifndef PLASK__SOLVER_VSLAB_MATRIX_H
18#define PLASK__SOLVER_VSLAB_MATRIX_H
19
20#include <cstring>
21#include <cmath>
22
23#include <plask/plask.hpp>
24#include "fortran.hpp"
25
26 namespace plask { namespace optical { namespace modal {
27
28 template <typename T> class MatrixDiagonal;
29
31 template <typename T>
32 class Matrix {
33 protected:
34 size_t r, c;
35
36 T* data_;
37 std::atomic<int>* gc;
38
39 void dec_ref() { // see http://www.boost.org/doc/libs/1_53_0/doc/html/atomic/usage_examples.html "Reference counting" for optimal memory access description
40 if (gc && gc->fetch_sub(1, std::memory_order_release) == 1) {
41 std::atomic_thread_fence(std::memory_order_acquire);
42 delete gc;
43 aligned_free(data_);
44 write_debug("freeing matrix {:d}x{:d} ({:.3f} MB) at {:p}", r, c, double(r*c*sizeof(T))/1048576., (void*)data_);
45 }
46 }
47
48 void inc_ref() {
49 if (gc) gc->fetch_add(1, std::memory_order_relaxed);
50 }
51
52 public:
53
54 typedef T DataType;
55
56 Matrix() : r(0), c(0), data_(nullptr), gc(nullptr) {}
57
58 Matrix(size_t m, size_t n) : r(m), c(n),
59 data_(aligned_malloc<T>(m*n)), gc(new std::atomic<int>(1)) {
60 write_debug("allocating matrix {:d}x{:d} ({:.3f} MB) at {:p}", r, c, double(r*c*sizeof(T))/1048576., (void*)data_);
61 }
62
63 Matrix(size_t m, size_t n, T val) : r(m), c(n),
64 data_(aligned_malloc<T>(m*n)), gc(new std::atomic<int>(1)) {
65 write_debug("allocating matrix {:d}x{:d} ({:.3f} MB) at {:p}", r, c, double(r*c*sizeof(T))/1048576., (void*)data_);
66 std::fill_n(data_, m*n, val);
67 }
68
69 Matrix(const Matrix<T>& M) : r(M.r), c(M.c), data_(M.data_), gc(M.gc) {
70 inc_ref();
71 }
72
73 Matrix<T>& operator=(const Matrix<T>& M) {
74 const_cast<Matrix<T>&>(M).inc_ref(); // must be called before dec_ref in case of self-asigment with *gc==1!
75 this->dec_ref();
76 r = M.r; c = M.c; data_ = M.data_; gc = M.gc;
77 return *this;
78 }
79
80 Matrix(const MatrixDiagonal<T>& M): r(M.size()), c(M.size()),
81 data_(aligned_malloc<T>(M.size()*M.size())), gc(new std::atomic<int>(1)) {
82 write_debug("allocating matrix {:d}x{:d} ({:.3f} MB) at {:p} (from diagonal)", r, c, double(r*c*sizeof(T))/1048576., (void*)data_);
83 std::fill_n(data_, r*c, 0);
84 for (int j = 0, n = 0; j < r; j++, n += c+1) data_[n] = M[j];
85 }
86
87 Matrix(size_t m, size_t n, T* existing_data) : r(m), c(n), data_(existing_data), gc(nullptr) {
88 // Create matrix using exiting data. This data is never garbage-collected
89 }
90
91 Matrix(size_t m, size_t n, std::initializer_list<T> data): r(m), c(n),
92 data_(aligned_malloc<T>(m*n)), gc(new std::atomic<int>(1)) {
93 std::copy(data.begin(), data.end(), data_);
94 }
95
96 ~Matrix() {
97 dec_ref();
98 }
99
100 void reset() {
101 dec_ref();
102 r = c = 0;
103 data_ = nullptr;
104 gc = nullptr;
105 }
106
107 void reset(size_t m, size_t n) {
108 dec_ref();
109 r = m; c = n;
110 data_ = aligned_malloc<T>(m*n);
111 gc = new std::atomic<int>(1);
112 write_debug("allocating matrix {:d}x{:d} ({:.3f} MB) at {:p}", r, c, double(r*c*sizeof(T))/1048576., (void*)data_);
113 }
114
115 void reset(size_t m, size_t n, T val) {
116 dec_ref();
117 r = m; c = n;
118 data_ = aligned_malloc<T>(m*n);
119 gc = new std::atomic<int>(1);
120 write_debug("allocating matrix {:d}x{:d} ({:.3f} MB) at {:p}", r, c, double(r*c*sizeof(T))/1048576., (void*)data_);
121 std::fill_n(data_, m*n, val);
122 }
123
124 void reset(size_t m, size_t n, T* existing_data) {
125 dec_ref();
126 r = m; c = n;
127 data_ = existing_data;
128 gc = nullptr;
129 }
130
131 inline const T* data() const { return data_; }
132 inline T* data() { return data_; }
133
134 inline const T& operator[](size_t i) const {
135 assert(i < r*c);
136 return data_[i];
137 }
138 inline T& operator[](size_t i) {
139 assert(i < r*c);
140 return data_[i];
141 }
142
143 inline const T& operator()(size_t m, size_t n) const {
144 assert(m < r);
145 assert(n < c);
146 return data_[n*r + m];
147 }
148 inline T& operator()(size_t m, size_t n) {
149 assert(m < r);
150 assert(n < c);
151 return data_[n*r + m];
152 }
153
154 inline size_t rows() const { return r; }
155 inline size_t cols() const { return c; }
156
157 Matrix<T> copy() const {
158 Matrix<T> copy_(r, c);
159 std::copy_n(data_, r*c, copy_.data());
160 return copy_;
161 }
162
163 void copyto(Matrix<T>& dst) {
164 assert(dst.rows() == rows() && dst.cols() == cols());
165 assert(dst.rows() == rows() && dst.cols() == cols());
166 std::copy(begin(), end(), dst.begin());
167 }
168
169 Matrix<T>& operator*=(T a) {
170 size_t size = r*c; for (size_t i = 0; i < size; i++) data_[i] *= a;
171 return *this;
172 }
173 Matrix<T>& operator/=(T a) {
174 size_t size = r*c; for (size_t i = 0; i < size; i++) data_[i] /= a;
175 return *this;
176 }
177
179 inline bool isnan() const {
180 const size_t n = r * c;
181 for (size_t i = 0; i < n; ++i)
182 if (std::isnan(real(data_[i])) || std::isnan(imag(data_[i]))) return true;
183 return false;
184 }
185
187 bool empty() const {
188 return !data_;
189 }
190
191 T* begin() const {
192 return data_;
193 }
194
195 T* end() const {
196 return data_ + r*c;
197 }
198};
199
200
202 template <typename T>
203 class MatrixDiagonal {
204 protected:
205 size_t siz;
206
207 T* data_; //< The data of the matrix
208 std::atomic<int>* gc; //< the reference count for the garbage collector
209
210 void dec_ref() { // see http://www.boost.org/doc/libs/1_53_0/doc/html/atomic/usage_examples.html "Reference counting" for optimal memory access description
211 if (gc && gc->fetch_sub(1, std::memory_order_release) == 1) {
212 std::atomic_thread_fence(std::memory_order_acquire);
213 delete gc;
214 aligned_free(data_);
215 write_debug("freeing diagonal matrix {0}x{0} ({1:.3f} MB) at {2}", siz, double(siz*sizeof(T))/1048576., (void*)data_);
216 }
217 }
218
219 void inc_ref() {
220 if (gc) gc->fetch_add(1, std::memory_order_relaxed);
221 }
222
223 public:
224
225 typedef T DataType;
226
227 MatrixDiagonal() : siz(0), gc(nullptr) {}
228
229 MatrixDiagonal(size_t n) : siz(n), data_(aligned_malloc<T>(n)), gc(new std::atomic<int>(1)) {
230 write_debug("allocating diagonal matrix {0}x{0} ({1:.3f} MB) at {2}", siz, double(siz*sizeof(T))/1048576., (void*)data_);
231 }
232
233 MatrixDiagonal(size_t n, T val) : siz(n), data_(aligned_malloc<T>(n)), gc(new std::atomic<int>(1)) {
234 write_debug("allocating and filling diagonal matrix {0}x{0} ({1:.3f} MB) at {2}", siz, double(siz*sizeof(T))/1048576., (void*)data_);
235 std::fill_n(data_, n, val);
236 }
237
238 MatrixDiagonal(const MatrixDiagonal<T>& M) : siz(M.siz), data_(M.data_), gc(M.gc) {
239 inc_ref();
240 }
241
242 MatrixDiagonal<T>& operator=(const MatrixDiagonal<T>& M) {
243 const_cast<MatrixDiagonal<T>&>(M).inc_ref(); // must be called before dec_ref in case of self-asigment with *gc==1!
244 this->dec_ref();
245 siz = M.siz; data_ = M.data_; gc = M.gc;
246 return *this;
247 }
248
249 ~MatrixDiagonal() {
250 dec_ref();
251 }
252
253 inline void reset() {
254 dec_ref();
255 siz = 0;
256 data_ = nullptr;
257 gc = nullptr;
258 }
259
260 void reset(size_t n) {
261 dec_ref();
262 siz = n;
263 data_ = aligned_malloc<T>(n);
264 gc = new std::atomic<int>(1);
265 write_debug("allocating diagonal matrix {0}x{0} ({1:.3f} MB) at {2}", siz, double(siz*sizeof(T))/1048576., (void*)data_);
266 }
267
268 void reset(size_t n, T val) {
269 dec_ref();
270 siz = n;
271 data_ = aligned_malloc<T>(n);
272 gc = new std::atomic<int>(1);
273 write_debug("allocating diagonal matrix {0}x{0} ({1:.3f} MB) at {2}", siz, double(siz*sizeof(T))/1048576., (void*)data_);
274 std::fill_n(data_, n, val);
275 }
276
277 inline const T* data() const { return data_; }
278 inline T* data() { return data_; }
279
280 inline const T& operator[](size_t n) const {
281 assert(n < siz);
282 return data_[n];
283 }
284 inline T& operator[](size_t n) {
285 assert(n < siz);
286 return data_[n];
287 }
288
289 inline const T& operator()(size_t m, size_t n) const {
290 assert(m < siz);
291 assert(n < siz);
292 return (m == n)? data_[n] : 0;
293 }
294 inline T& operator()(size_t m, size_t n) {
295 assert(m < siz);
296 assert(n < siz);
297 assert(m == n);
298 return data_[n];
299 }
300
301 inline size_t size() const { return siz; }
302
303 MatrixDiagonal<T> copy() const {
304 MatrixDiagonal<T> C(siz);
305 std::copy_n(data_, siz, C.data());
306 return C;
307 }
308
309 MatrixDiagonal<T>& operator*=(T a) { for (int i = 0; i < siz; i++) data_[i] *= a; return *this; }
310 MatrixDiagonal<T>& operator/=(T a) { for (int i = 0; i < siz; i++) data_[i] /= a; return *this; }
311
313 inline bool isnan() const {
314 for (size_t i = 0; i != siz; ++i)
315 if (std::isnan(real(data_[i])) || std::isnan(imag(data_[i]))) return true;
316 return false;
317 }
318
320 bool empty() const {
321 return !data_;
322 }
323
324 T* begin() const {
325 return data_;
326 }
327
328 T* end() const {
329 return data_ + siz;
330 }
331};
332
333//**************************************************************************
334// Rectangular matrix of real and complex numbers
335 typedef Matrix<double> dmatrix;
336 typedef Matrix<dcomplex> cmatrix;
337
338// Column vector and diagonal matrix
339 typedef DataVector<dcomplex> cvector;
340 typedef DataVector<double> dvector;
341 typedef DataVector<const dcomplex> const_cvector;
342 typedef MatrixDiagonal<dcomplex> cdiagonal;
343
344//**************************************************************************
346 template <typename T>
347 inline Matrix<T> operator+(const Matrix<T>& A, const Matrix<T>& B) {
348 // if (A.rows() != B.rows() || A.cols() != B.cols()) throw ComputationError("operator*<cmatrix,cmatrix>", "cannot add: dimensions does not match");
349 assert(A.rows() == B.rows() && A.cols() == B.cols());
350 Matrix<T> C(A.rows(), A.cols());
351 T* c = C.data();
352 for (const T *a = A.data(), *b = B.data(), *end = A.data() + A.rows()*A.cols(); a != end; ++a, ++b, ++c)
353 *c = *a + *b;
354 return C;
355}
356
358 template <typename T>
359 inline Matrix<T> operator-(const Matrix<T>& A, const Matrix<T>& B) {
360 // if (A.rows() != B.rows() || A.cols() != B.cols()) throw ComputationError("operator*<cmatrix,cmatrix>", "cannot add: dimensions does not match");
361 assert(A.rows() == B.rows() && A.cols() == B.cols());
362 Matrix<T> C(A.rows(), A.cols());
363 T* c = C.data();
364 for (const T *a = A.data(), *b = B.data(), *end = A.data() + A.rows()*A.cols(); a != end; ++a, ++b, ++c)
365 *c = *a - *b;
366 return C;
367}
368
370 template <typename T>
371 inline Matrix<T> operator-(const Matrix<T>& A) {
372 Matrix<T> C(A.rows(), A.cols());
373 for (T *a = A.data(), *c = C.data(), *end = A.data() + A.rows()*A.cols(); a != end; ++a, ++c)
374 *c = - *a;
375 return C;
376}
377
379 inline cmatrix operator*(const cmatrix& A, const cmatrix& B) {
380 // if (A.cols() != B.rows()) throw ComputationError("operator*<cmatrix,cmatrix>", "cannot multiply: A.cols != B.rows");
381 assert(A.cols() == B.rows());
382 cmatrix C(A.rows(), B.cols());
383 zgemm('n', 'n', int(A.rows()), int(B.cols()), int(A.cols()), 1., A.data(), int(A.rows()), B.data(), int(B.rows()), 0., C.data(), int(C.rows()));
384 return C;
385}
386
388 inline cvector operator*(const cmatrix& A, const cvector& v) {
389 size_t n = A.cols(), m = A.rows();
390 // if (n != v.size()) throw ComputationError("mult_matrix_by_vector", "A.cols != v.size");
391 assert(n == v.size());
392 cvector dst(m);
393 zgemv('n', int(m), int(n), 1., A.data(), int(m), v.data(), 1, 0., dst.data(), 1);
394 return dst;
395}
396
398 template <typename T>
399 inline cmatrix operator*(const Matrix<T>& A, const MatrixDiagonal<T>& B) {
400 // if (A.cols() != B.size()) throw ComputationError("operator*<cmatrix,cdiagonal>", "cannot multiply: A.cols != B.size");
401 assert(A.cols() == B.size());
402 cmatrix C(A.rows(), B.size());
403 size_t n = 0;
404 const size_t r = A.rows(), c = A.cols();
405 for (size_t j = 0; j < c; j++)
406 {
407 T b = B[j];
408 for (size_t i = 0; i < r; i++, n++)
409 C[n] = A[n] * b;
410 }
411 return C;
412}
413
415 template <typename T>
416 inline cmatrix operator*(const MatrixDiagonal<T>& A, const Matrix<T>& B) {
417 // if (A.size() != B.rows()) throw ComputationError("operator*<cdiagonal,cmatrix>", "cannot multiply: A.size != B.rows");
418 assert(A.size() == B.rows());
419 cmatrix C(A.size(), B.cols());
420 size_t n = 0;
421 const size_t r = B.rows(), c = B.cols();
422 for (size_t j = 0; j < c; j++)
423 for (size_t i = 0; i < r; i++, n++)
424 C[n] = A[i] * B[n];
425 return C;
426}
427
429 template <typename T>
430 inline void mult_matrix_by_diagonal(Matrix<T>& A, const MatrixDiagonal<T>& B) {
431 // if (A.cols() != B.size()) throw ComputationError("mult_matrix_by_diagonal", "cannot multiply: A.cols != B.size");
432 assert(A.cols() == B.size());
433 size_t n = 0;
434 const size_t r = A.rows(), c = A.cols();
435 for (size_t j = 0; j < c; j++) {
436 T b = B[j];
437 for (size_t i = 0; i < r; i++, n++)
438 A[n] *= b;
439 }
440}
441
443 template <typename T>
444 inline void mult_diagonal_by_matrix(const MatrixDiagonal<T>& A, Matrix<T>& B) {
445 // if (A.size() != B.rows()) throw ComputationError("mult_diagonal_by_matrix", "cannot multiply: A.size != B.rows");
446 assert(A.size() == B.rows());
447 size_t n = 0;
448 const size_t r = B.rows(), c = B.cols();
449 for (size_t j = 0; j < c; j++)
450 for (size_t i = 0; i < r; i++, n++)
451 B[n] *= A[i];
452}
453
454
455// BLAS wrappers for multiplications without allocating additional storage
456 inline void mult_matrix_by_vector(const cmatrix& A, const const_cvector& v, cvector& dst) {
457 const size_t m = A.rows(),
458 n = A.cols();
459 // if (n != v.size()) throw ComputationError("mult_matrix_by_vector", "A.cols != v.size");
460 // if (m != dst.size()) throw ComputationError("mult_matrix_by_vector", "A.rows != dst.size");
461 assert(n == v.size());
462 assert(m == dst.size());
463 zgemv('n', int(m), int(n), 1., A.data(), int(m), v.data(), 1, 0., dst.data(), 1);
464}
465
466 inline void mult_matrix_by_matrix(const cmatrix& A, const cmatrix& B, cmatrix& dst) {
467 const size_t k = A.cols(),
468 m = A.rows(),
469 n = B.cols();
470 // if (k != B.rows()) throw ComputationError("mult_matrix_by_matrix", "cannot multiply: A.cols != B.rows");
471 // if (m != dst.rows()) throw ComputationError("mult_matrix_by_matrix", "A.rows != dst.rows");
472 // if (n != dst.cols()) throw ComputationError("mult_matrix_by_matrix", "B.cols != dst.cols");
473 assert(k == B.rows());
474 assert(m == dst.rows());
475 assert(n == dst.cols());
476 zgemm('n', 'n', int(m), int(n), int(k), 1., A.data(), int(m), B.data(), int(k), 0., dst.data(), int(m));
477}
478
479 inline void add_mult_matrix_by_vector(const cmatrix& A, const cvector& v, cvector& dst) {
480 const size_t m = A.rows(),
481 n = A.cols();
482 // if (n != v.size()) throw ComputationError("add_mult_matrix_by_vector", "A.cols != v.size");
483 // if (m != dst.size()) throw ComputationError("add_mult_matrix_by_vector", "A.rows != dst.size");
484 assert(n == v.size());
485 assert(m == dst.size());
486 zgemv('n', int(m), int(n), 1., A.data(), int(m), v.data(), 1, 1., dst.data(), 1);
487}
488
489 inline void add_mult_matrix_by_matrix(const cmatrix& A, const cmatrix& B, cmatrix& dst) {
490 const size_t k = A.cols(),
491 m = A.rows(),
492 n = B.cols();
493 // if (k != B.rows()) throw ComputationError("add_mult_matrix_by_matrix", "cannot multiply: A.cols != B.rows");
494 // if (m != dst.rows()) throw ComputationError("add_mult_matrix_by_matrix", "A.rows != dst.rows");
495 // if (n != dst.cols()) throw ComputationError("add_mult_matrix_by_matrix", "B.cols != dst.cols");
496 assert(k == B.rows());
497 assert(m == dst.rows());
498 assert(n == dst.cols());
499 zgemm('n', 'n', int(m), int(n), int(k), 1., A.data(), int(m), B.data(), int(k), 1., dst.data(), int(m));
500}
501
502
503// Some LAPACK wrappers
504 cmatrix invmult(cmatrix& A, cmatrix& B);
505 cvector invmult(cmatrix& A, cvector& B);
506 cmatrix inv(cmatrix& A);
507dcomplex det(cmatrix& A);
508 int eigenv(cmatrix& A, cdiagonal& vals, cmatrix* rightv=NULL, cmatrix* leftv=NULL);
509
510// r-value wrappers
511 inline cmatrix invmult(cmatrix&& A, cmatrix& B) { return invmult(A, B); }
512 inline cvector invmult(cmatrix&& A, cvector& B) { return invmult(A, B); };
513 inline cmatrix inv(cmatrix&& A) { return inv(A); }
514 inline dcomplex det(cmatrix&& A) { return det(A); }
515 inline int eigenv(cmatrix&& A, cdiagonal& vals, cmatrix* rightv=NULL, cmatrix* leftv=NULL) {
516 return eigenv(A, vals, rightv, leftv);
517}
518
519
521 template <typename T>
522 void zero_matrix(Matrix<T>& A, size_t m, size_t n) {
523 if (A.rows() != m || A.cols() != n) A.reset(m, n, 0.);
524 else std::fill_n(A.data(), m*n, T(0.));
525}
526
528 template <typename T>
529 void zero_matrix(Matrix<T>& A) {
530 std::fill_n(A.data(), A.rows()*A.cols(), T(0.));
531}
532
534 template <typename T>
535 void make_unit_matrix(Matrix<T>& A, size_t n) {
536 zero_matrix(A, n, n);
537 constexpr T one(1.);
538 for (int i = 0; i < n; ++i) A(i,i) = one;
539}
540
542 template <typename T>
543 void make_unit_matrix(Matrix<T>& A) {
544 assert(A.rows() == A.cols());
545 zero_matrix(A);
546 constexpr T one(1.);
547 for (int i = 0; i < A.rows(); ++i) A(i,i) = one;
548}
549
550}}}
551#endif // PLASK__SOLVER_VSLAB_MATRIX_H