PLaSK library

Loading...
Searching...
No Matches
diagonalizer.cpp
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#include "diagonalizer.hpp"
15#include "expansion.hpp"
16
17#include <plask/plask.hpp>
18
19#include <algorithm>
20#include <cstring>
21
22#ifdef OPENMP_FOUND
23# include <omp.h>
24#endif
25
26 namespace plask { namespace optical { namespace modal {
27
28 Diagonalizer::Diagonalizer(Expansion* src) :
29 src(src), diagonalized(src->solver->lcount, false), lcount(src->solver->lcount) {}
30
31 Diagonalizer::~Diagonalizer() {}
32
33
34 SimpleDiagonalizer::SimpleDiagonalizer(Expansion* g) :
35 Diagonalizer(g), gamma(lcount), Te(lcount), Th(lcount), Te1(lcount), Th1(lcount)
36{
37 const std::size_t N = src->matrixSize(); // Size of each matrix
38
39 for (std::size_t i = 0; i < lcount; i++) {
40 gamma[i].reset(N);
41 Th[i].reset(N, N);
42 Th1[i].reset(N, N);
43 Te[i].reset(N, N);
44 Te1[i].reset(N, N);
45 }
46}
47
48
49std::size_t SimpleDiagonalizer::matrixSize() const
50{
51 return src->matrixSize();
52}
53
54 void SimpleDiagonalizer::initDiagonalization()
55{
56 for (std::size_t layer = 0; layer < lcount; layer++)
57 diagonalized[layer] = false;
58}
59
60
61 bool SimpleDiagonalizer::diagonalizeLayer(size_t layer)
62{
63 if (diagonalized[layer]) return false;
64
65 const size_t N = src->matrixSize(); // Size of each matrix
66 const size_t NN = N*N;
67 cdiagonal& gam = gamma[layer];
68
69 #ifdef OPENMP_FOUND
70 writelog(LOG_DEBUG, "{}: Diagonalizing matrix for layer {:d}/{:d} in thread {:d}",
71 src->solver->getId(), layer, lcount, omp_get_thread_num());
72 #else
73 writelog(LOG_DEBUG, "{}: Diagonalizing matrix for layer {:d}/{:d}", src->solver->getId(), layer, lcount);
74 #endif
75
76 // First find necessary matrices
77 cmatrix RE = Th1[layer], RH = Th[layer];
78
79 src->getMatrices(layer, RE, RH);
80
81 // Ugly hack to avoid singularities
82 for (std::size_t i = 0; i != N; ++i) {
83 if (RE(i,i) == 0.) RE(i,i) = SMALL;
84 if (RH(i,i) == 0.) RH(i,i) = SMALL;
85 }
86
87 // std::cerr << "PLaSK\nRE:\n";
88 // for (unsigned r = 0; r != N; ++r) {
89 // for (unsigned c = 0; c != N; ++c)
90 // std::cerr << format("{:7.1f} ", real(RE(r,c)));
91 // std::cerr << "\n";
92 // }
93 // std::cerr << "RH:\n";
94 // for (unsigned r = 0; r != N; ++r) {
95 // for (unsigned c = 0; c != N; ++c)
96 // std::cerr << format("{:7.1f} ", real(RH(r,c)));
97 // std::cerr << "\n";
98 // }
99 assert(!RE.isnan());
100 assert(!RH.isnan());
101
102 TempMatrix temp = src->getTempMatrix();
103 cmatrix QE(temp);
104
105 if (src->diagonalQE(layer)) {
106
107 // We are lucky - the QH matrix is diagonal so we can make it fast and easy
108
109 // So we compute the diagonal elements of QH = RE*RH
110 for (std::size_t ie = 0, ih = 0; ie < N; ie++, ih += N) {
111 gam[ie] = 0;
112 for (std::size_t jh = 0, je = 0; jh < N; jh++, je += N)
113 gam[ie] += RH[ie+je] * RE[ih+jh];
114 }
115
116 // Make Gamma of Gamma^2
117 sqrtGamma(gam);
118
119 src->getDiagonalEigenvectors(Te[layer], Te1[layer], RE, gam);
120
121 } else {
122 // We have to make the proper diagonalization
123 // TODO: rewrite it to more low-level and more optimized computations
124 mult_matrix_by_matrix(RH, RE, QE); // QE = RH * RE
125
126 // std::cerr << "PLaSK\nQE:\n";
127 // for (unsigned r = 0; r != N; ++r) {
128 // for (unsigned c = 0; c != N; ++c)
129 // std::cerr << format("{:7.1f} ", real(QE(r,c)));
130 // std::cerr << "\n";
131 // }
132
133 // This is probably expensive but necessary check to avoid hangs
134 if (QE.isnan()) throw ComputationError(src->solver->getId(), "simpleDiagonalizer: NaN in Q matrix");
135
136 // Here we make the actual diagonalization, i.e. compute the eigenvalues and eigenvectors of QE
137 int info;
138 if (N < 2) {
139 dcomplex lwork[4];
140 double rwork[2];
141 zgeev('N', 'V', int(N), QE.data(), int(N), gam.data(), nullptr, int(N), Te[layer].data(), int(N),
142 lwork, 2, rwork, info);
143 } else {
144 // We use Th as work and Te1 as rwork (as N >= 2, their sizes are ok)
145 zgeev('N', 'V', int(N), QE.data(), int(N), gam.data(), nullptr, int(N), Te[layer].data(), int(N),
146 Th[layer].data(), int(NN), reinterpret_cast<double*>(Te1[layer].data()), info);
147 }
148 if (info != 0) throw ComputationError(src->solver->getId(), "simpleDiagonalizer: Could not compute {0}-th eignevalue of QE", info);
149
150 // Find the inverse of Te in the classical way (maybe to be optimized in future)
151 // TODO: eigenvectors should be built by hand based on Schur vectors
152 std::copy_n(Te[layer].data(), NN, Th[layer].data());
153 make_unit_matrix(Te1[layer]);
154 invmult(Th[layer], Te1[layer]);
155
156 // Make Gamma of Gamma^2
157 sqrtGamma(gam);
158 }
159 assert(!Te[layer].isnan());
160
161 // So now there is the time to find TH = Re * Te * Gamma^(-1)
162 mult_matrix_by_matrix(RE, Te[layer], Th[layer]);
163 dcomplex* th = Th[layer].data();
164 for (std::size_t j = 0; j < N; j++) {
165 dcomplex g = 1. / gam[j];
166 for (std::size_t i = 0; i < N; i++) *(th+i) *= g;
167 th += N;
168 }
169 assert(!Th[layer].isnan());
170
171 // Compute the Th1[layer] = Gamma[layer] * Te1[layer] * inv(RE)
172 // we use the LU factorization of the RE matrix for this purpose and then solve Th1^T = inv(RE)^T * Te1^T * Gamma^T
173 // the QE array is used as a temporary storage
174 for (std::size_t i = 0; i < N; i++)
175 for (std::size_t j = 0; j < N; j++)
176 QE(i,j) = Te1[layer](j,i);
177 // LU factorization of RE
178 int ierr;
179 std::unique_ptr<int[]> ipiv(new int[N]);
180 zgetrf(int(N), int(N), RE.data(), int(N), ipiv.get(), ierr);
181 if (ierr != 0) throw ComputationError(src->solver->getId(), "simpleDiagonalizer: RE matrix singular");
182 // the QE will contain inv(RE)^T * Te1^T
183 zgetrs('t', int(N), int(N), RE.data(), int(N), ipiv.get(), QE.data(), int(N), ierr);
184 if (ierr != 0) throw ComputationError(src->solver->getId(), "simpleDiagonalizer: Could not compute inv(RE)");
185 // compute QE^T and store it in Th1
186 for (std::size_t j = 0; j < N; j++) {
187 dcomplex g = gam[j];
188 for (std::size_t i = 0; i < N; i++)
189 Th1[layer](j,i) = QE(i,j) * g;
190 }
191 assert(!Th1[layer].isnan());
192
193 // Mark that layer has been diagonalized
194 diagonalized[layer] = true;
195
196 return true;
197}
198
199}}} // namespace plask::optical::modal