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