1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
// Test ../include/Spectra/LinAlg/UpperHessenbergEigen.h and
// ../include/Spectra/LinAlg/TridiagEigen.h
#include <chrono>
#include <iostream>
#include <complex>
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
#include <Spectra/LinAlg/UpperHessenbergEigen.h>
#include <Spectra/LinAlg/TridiagEigen.h>
#include "catch.hpp"
using namespace Spectra;
// https://stackoverflow.com/a/34781413
using Clock = std::chrono::high_resolution_clock;
using Duration = std::chrono::duration<double, std::milli>;
using TimePoint = std::chrono::time_point<Clock, Duration>;
using Matrix = Eigen::MatrixXd;
using Vector = Eigen::VectorXd;
using Complex = std::complex<double>;
using ComplexMatrix = Eigen::MatrixXcd;
using ComplexVector = Eigen::VectorXcd;
TEST_CASE("Eigen decomposition of real upper Hessenberg matrix", "[Eigen]")
{
std::srand(123);
const int n = 100;
Matrix M = Matrix::Random(n, n);
// H is upper Hessenberg
Matrix H = M.triangularView<Eigen::Upper>();
H.diagonal(-1) = M.diagonal(-1);
UpperHessenbergEigen<double> decomp(H);
ComplexVector evals = decomp.eigenvalues();
ComplexMatrix evecs = decomp.eigenvectors();
// Test accuracy
ComplexMatrix err = H * evecs - evecs * evals.asDiagonal();
INFO("||HU - UD||_inf = " << err.cwiseAbs().maxCoeff());
REQUIRE(err.cwiseAbs().maxCoeff() == Approx(0.0).margin(1e-12));
TimePoint t1, t2;
t1 = Clock::now();
for (int i = 0; i < 100; i++)
{
UpperHessenbergEigen<double> decomp(H);
ComplexVector evals = decomp.eigenvalues();
ComplexMatrix evecs = decomp.eigenvectors();
}
t2 = Clock::now();
std::cout << "Elapsed time for UpperHessenbergEigen: "
<< (t2 - t1).count() << " ms\n";
t1 = Clock::now();
for (int i = 0; i < 100; i++)
{
Eigen::EigenSolver<Matrix> decomp(H);
ComplexVector evals = decomp.eigenvalues();
ComplexMatrix evecs = decomp.eigenvectors();
}
t2 = Clock::now();
std::cout << "Elapsed time for Eigen::EigenSolver: "
<< (t2 - t1).count() << " ms\n";
}
TEST_CASE("Eigen decomposition of real symmetric tridiagonal matrix", "[Eigen]")
{
std::srand(123);
const int n = 100;
Matrix M = Matrix::Random(n, n);
// H is symmetric tridiagonal
Matrix H = Matrix::Zero(n, n);
H.diagonal() = M.diagonal();
H.diagonal(-1) = M.diagonal(-1);
H.diagonal(1) = M.diagonal(-1);
TridiagEigen<double> decomp(H);
Vector evals = decomp.eigenvalues();
Matrix evecs = decomp.eigenvectors();
// Test accuracy
Matrix err = H * evecs - evecs * evals.asDiagonal();
INFO("||HU - UD||_inf = " << err.cwiseAbs().maxCoeff());
REQUIRE(err.cwiseAbs().maxCoeff() == Approx(0.0).margin(1e-12));
TimePoint t1, t2;
t1 = Clock::now();
for (int i = 0; i < 100; i++)
{
TridiagEigen<double> decomp(H);
Vector evals = decomp.eigenvalues();
Matrix evecs = decomp.eigenvectors();
}
t2 = Clock::now();
std::cout << "Elapsed time for TridiagEigen: "
<< (t2 - t1).count() << " ms\n";
t1 = Clock::now();
for (int i = 0; i < 100; i++)
{
Eigen::SelfAdjointEigenSolver<Matrix> decomp(H);
Vector evals = decomp.eigenvalues();
Matrix evecs = decomp.eigenvectors();
}
t2 = Clock::now();
std::cout << "Elapsed time for Eigen::SelfAdjointEigenSolver: "
<< (t2 - t1).count() << " ms\n";
}
TEST_CASE("Eigen decomposition of complex upper Hessenberg matrix", "[Eigen]")
{
std::srand(123);
const int n = 100;
ComplexMatrix M = ComplexMatrix::Random(n, n);
// H is upper Hessenberg
ComplexMatrix H = M.triangularView<Eigen::Upper>();
H.diagonal(-1) = M.diagonal(-1);
UpperHessenbergEigen<Complex> decomp(H);
ComplexVector evals = decomp.eigenvalues();
ComplexMatrix evecs = decomp.eigenvectors();
// Test accuracy
ComplexMatrix err = H * evecs - evecs * evals.asDiagonal();
INFO("||HU - UD||_inf = " << err.cwiseAbs().maxCoeff());
REQUIRE(err.cwiseAbs().maxCoeff() == Approx(0.0).margin(1e-12));
TimePoint t1, t2;
t1 = Clock::now();
for (int i = 0; i < 100; i++)
{
UpperHessenbergEigen<Complex> decomp(H);
ComplexVector evals = decomp.eigenvalues();
ComplexMatrix evecs = decomp.eigenvectors();
}
t2 = Clock::now();
std::cout << "Elapsed time for UpperHessenbergEigen: "
<< (t2 - t1).count() << " ms\n";
t1 = Clock::now();
for (int i = 0; i < 100; i++)
{
Eigen::ComplexEigenSolver<ComplexMatrix> decomp(H);
ComplexVector evals = decomp.eigenvalues();
ComplexMatrix evecs = decomp.eigenvectors();
}
t2 = Clock::now();
std::cout << "Elapsed time for Eigen::EigenSolver: "
<< (t2 - t1).count() << " ms\n";
}
|
