PLaSK library

Loading...
Searching...
No Matches
gauss_matrix.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 */
14#ifndef PLASK__SOLVER__OPTICAL_EFFECTIVE_GAUSS_MATRIX_H
15#define PLASK__SOLVER__OPTICAL_EFFECTIVE_GAUSS_MATRIX_H
16
17#include <cstddef>
18#include <plask/plask.hpp>
19 using plask::dcomplex;
20
21// BLAS routine to multiply matrix by vector
22#define zgbmv F77_GLOBAL(zgbmv,ZGBMV)
23F77SUB zgbmv(const char& trans, const int& m, const int& n, const int& kl, const int& ku, const dcomplex& alpha, dcomplex* a, const int& lda,
24 const dcomplex* x, int incx, const dcomplex& beta, dcomplex* y, int incy);
25
26
27// LAPACK routines to solve set of linear equations
28#define zgbtrf F77_GLOBAL(zgbtrf,ZGBTRF)
29F77SUB zgbtrf(const int& m, const int& n, const int& kl, const int& ku, dcomplex* ab, const int& ldab, int* ipiv, int& info);
30
31#define zgbtrs F77_GLOBAL(zgbtrs,ZGBTRS)
32F77SUB zgbtrs(const char& trans, const int& n, const int& kl, const int& ku, const int& nrhs, dcomplex* ab, const int& ldab, int* ipiv, dcomplex* b, const int& ldb, int& info);
33
34 namespace plask { namespace optical { namespace effective {
35
36 constexpr const int LD = 7;
37
42 struct ZgbMatrix {
43
44 const size_t size;
45 dcomplex* data;
46
51 ZgbMatrix(size_t rank): size(rank), data(aligned_malloc<dcomplex>(LD*rank)) {}
52
53 ZgbMatrix(const ZgbMatrix&) = delete; // this object is non-copyable
54
55 ~ZgbMatrix() { aligned_free(data); }
56
62 size_t index(size_t r, size_t c) {
63 assert(r < size && c < size);
64 assert(abs(int(c)-int(r)) < 3);
65 // AB(kl+ku+1+i-j,j) = A(i,j)
66 return (LD-1)*c + r + 4;
67 }
68
74 dcomplex& operator()(size_t r, size_t c) {
75 return data[index(r,c)];
76 }
77
79 void clear() {
80 std::fill_n(data, LD*size, 0.);
81 }
82
88 void mult(const DataVector<const dcomplex>& vector, DataVector<dcomplex>& result) {
89 zgbmv('N', int(size), int(size), 2, 2, 1., data, LD, vector.data(), 1, 0., result.data(), 1);
90 }
91
97 void addmult(const DataVector<const dcomplex>& vector, DataVector<dcomplex>& result) {
98 zgbmv('N', int(size), int(size), 2, 2, 1., data, LD, vector.data(), 1, 1., result.data(), 1);
99 }
100
102 dcomplex determinant() {
103 int info = 0;
104 aligned_unique_ptr<int> upiv(aligned_malloc<int>(size));
105 int* ipiv = upiv.get();
106 zgbtrf(int(size), int(size), 2, 2, data, LD, ipiv, info);
107 assert(info >= 0);
108
109 dcomplex det = 1.;
110 for (std::size_t i = 0; i < size; ++i) {
111 det *= data[LD*i + 4];
112 if (ipiv[i] != int(i+1)) det = -det;
113 }
114 return det;
115 }
116};
117
118}}} // # namespace plask::gain::freecarrier
119
120#endif // PLASK__SOLVER__OPTICAL_EFFECTIVE_GAUSS_MATRIX_H
121