PLaSK library

Loading...
Searching...
No Matches
invmult_benchmark.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 <plask/plask.hpp>
15#include "../matrices.hpp"
16#include "../fourier/toeplitz.hpp"
17
18 using namespace plask;
19 using namespace plask::optical::modal;
20
21
22 const int n = 1000;
23
24
25 const int N = 2*n + 1;
26 const int nN = 4*n + 1;
27 const double b = 2*PI;
28
29 const dcomplex e(3., 0.1);
30
31 int main() {
32 DataVector<dcomplex> eps(nN);
33
34 for (int k = -N+1; k < N; ++k) {
35 size_t j = (k>=0)? k : k + nN;
36 dcomplex ff = (j)? (dcomplex(0., 0.5/PI/k) * (exp(dcomplex(0., -b*k*0.25)) - exp(dcomplex(0., +b*k*0.15)))) : 0.5;
37 eps[j] = e * ff;
38 // eps[j] = dcomplex(k, 0.);
39 }
40 eps[0] += 1.;
41
42// for(size_t j = 0; j < nN; ++j)
43// std::cerr << str(eps[j], "{:9.6f}{:+09.6f}j") << " ";
44// std::cerr << "\n\n";
45
46 cmatrix T(N, N);
47
48// dcomplex* p = T.data();
49// for (size_t j = 0; j < N; ++j, p += N+1) {
50// std::copy_n(eps.data(), N-j, p); // below diagonal
51// std::copy_n(eps.data()+nN-j, j, T.data()+N*j); // above diagonal
52// }
53 for (int i = -n; i <= n; ++i) {
54 for (int j = -n; j <= n; ++j) {
55 int ij = i-j; if (ij < 0) ij += nN;
56 T((i>=0)?i:i+N, (j>=0)?j:j+N) = eps[ij];
57 // T(i+n, j+n) = eps[ij];
58 }
59 }
60
61 cmatrix X(N, N, 0.);
62 for (size_t i = 0; i < N; ++i) {
63 X(i, i) = 1.;
64 }
65
66// std::cerr << "T = matrix([\n";
67// for (size_t i = 0; i < N; ++i) {
68// std::cerr << " [ ";
69// for(size_t j = 0; j < N; ++j) {
70// std::cerr << str(T(i,j), "{:9.6f}{:+09.6f}j") << ", ";
71// }
72// std::cerr << "],\n";
73// }
74// std::cerr << "])\n\n";
75
76 // Invert matrices with LAPACK
77 std::unique_ptr<int[]> ipiv(new int[N]);
78 int info;
79 zgesv(N, N, T.data(), N, ipiv.get(), X.data(), N, info);
80 if (info > 0) throw ComputationError("invmult", "Toeplitz matrix singular");
81
82 // ToeplitzLevinson(eps, X);
83
84 // std::cerr << "X = matrix([\n";
85 // for (size_t i = 0; i < N; ++i) {
86 // std::cerr << " [ ";
87 // for(size_t j = 0; j < N; ++j) {
88 // std::cerr << str(X(i,j), "{:9.6f}{:+09.6f}j") << ", ";
89 // }
90 // std::cerr << "],\n";
91 // }
92 // std::cerr << "])\n";
93
94}