PLaSK library

Loading...
Searching...
No Matches
diagonalizer.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_OPTICAL_MODAL_DIAGONALIZER_H
18#define PLASK__SOLVER_OPTICAL_MODAL_DIAGONALIZER_H
19
20#include <utility>
21
22#ifdef OPENMP_FOUND
23# include <omp.h>
24#endif
25
26#include <plask/plask.hpp>
27
28#include "matrices.hpp"
29#include "expansion.hpp"
30
31 namespace plask { namespace optical { namespace modal {
32
42 class Diagonalizer
43{
44 protected:
45 Expansion* src;
46 std::vector<bool> diagonalized;
47
48 public:
49 const std::size_t lcount;
50
51 Diagonalizer(Expansion* src);
52
53 virtual ~Diagonalizer();
54
56 virtual std::size_t matrixSize() const = 0;
57
59 inline const Expansion* source() const { return src; }
60
62 inline Expansion* source() { return src; }
63
65 virtual void initDiagonalization() = 0;
66
69 virtual bool diagonalizeLayer(size_t layer) = 0;
70
72 bool isDiagonalized(size_t layer) {
73 return diagonalized[layer];
74 }
75
77 virtual const cdiagonal& Gamma(size_t layer) const = 0;
78
80 virtual const cmatrix& TE(size_t layer) const = 0;
81
83 virtual const cmatrix& TH(size_t layer) const = 0;
84
86 virtual const cmatrix& invTE(size_t layer) const = 0;
87
89 virtual const cmatrix& invTH(size_t layer) const = 0;
90
91 // Diagonalization function to compute the smallest eigenvalue of the provided matrix
92// virtual dcomplex smallest_eigevalue(const cmatrix& M) = 0;
93};
94
95
101 class SimpleDiagonalizer : public Diagonalizer
102{
103 protected:
104 std::vector<cdiagonal> gamma;
105 std::vector<cmatrix> Te, Th;
106 std::vector<cmatrix> Te1, Th1;
107
110 void sqrtGamma(cdiagonal& gam) {
111 const size_t N = src->matrixSize();
112 for (std::size_t j = 0; j < N; j++) {
113 dcomplex g = sqrt(gam[j]);
114 if (g == 0.) g = SMALL; // Ugly hack to avoid singularity!
115 if (real(g) < -SMALL) g = -g;
116 if (imag(g) > SMALL) g = -g;
117 gam[j] = g;
118 }
119 }
120
121 public:
122 SimpleDiagonalizer(Expansion* g);
123 ~SimpleDiagonalizer() {}
124
125 std::size_t matrixSize() const override;
126
127 void initDiagonalization() override;
128
129 bool diagonalizeLayer(size_t layer) override;
130
131 // Functions returning references to calculated matrices
132 const cdiagonal& Gamma(size_t layer) const override { return gamma[layer]; }
133 const cmatrix& TE(size_t layer) const override { return Te[layer]; }
134 const cmatrix& TH(size_t layer) const override { return Th[layer]; }
135 const cmatrix& invTE(size_t layer) const override { return Te1[layer]; }
136 const cmatrix& invTH(size_t layer) const override { return Th1[layer]; }
137};
138
139}}} // namespace plask::optical::modal
140#endif // PLASK__SOLVER_OPTICAL_MODAL_DIAGONALIZER_H