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
|
// Example reported in Issue #144
// https://github.com/yixuan/spectra/issues/144
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
#include <Spectra/SymEigsSolver.h>
#include <Spectra/SymEigsShiftSolver.h>
#include <Spectra/MatOp/DenseSymMatProd.h>
#include <Spectra/MatOp/DenseSymShiftSolve.h>
#include "catch.hpp"
using namespace Spectra;
using Matrix = Eigen::MatrixXd;
using Vector = Eigen::VectorXd;
Matrix construct_cycle_laplacian(int n)
{
// Initialize the Laplacian matrix
Matrix L = Matrix::Zero(n, n);
// Add the matrix entries, iterating over the rows
for (int i = 0; i < n; i++)
{
L(i, i) = 1;
L(i, (i + n - 1) % n) = -0.5;
L(i, (i + 1) % n) = -0.5;
}
return L;
}
void run_test(int n, int k, int m)
{
const Matrix M = construct_cycle_laplacian(n);
// True eigenvalues
Eigen::SelfAdjointEigenSolver<Matrix> es(M);
Vector true_evals = es.eigenvalues();
// Largest eigenvalues
DenseSymMatProd<double> op(M);
SymEigsSolver<DenseSymMatProd<double>> eigs(op, k, m);
eigs.init();
int nconv = eigs.compute(SortRule::LargestMagn, 1000, 1e-15, SortRule::SmallestAlge);
int niter = eigs.num_iterations();
int nops = eigs.num_operations();
INFO("nconv = " << nconv);
INFO("niter = " << niter);
INFO("nops = " << nops);
REQUIRE(eigs.info() == CompInfo::Successful);
Vector evals = eigs.eigenvalues();
Matrix evecs = eigs.eigenvectors();
Matrix resid = M.selfadjointView<Eigen::Lower>() * evecs - evecs * evals.asDiagonal();
double err = resid.array().abs().maxCoeff();
INFO("||AU - UD||_inf = " << err);
REQUIRE(err == Approx(0.0).margin(1e-9));
INFO("True eigenvalues =\n " << true_evals);
INFO("Estimated =\n " << evals);
double diff = (true_evals.tail(k) - evals).array().abs().maxCoeff();
INFO("diff = " << diff);
REQUIRE(diff == Approx(0.0).margin(1e-9));
// Smallest eigenvalues
DenseSymShiftSolve<double> op2(M);
SymEigsShiftSolver<DenseSymShiftSolve<double>> eigs2(op2, k, m, -1e-6);
eigs2.init();
nconv = eigs2.compute(SortRule::LargestMagn, 1000, 1e-15, SortRule::SmallestAlge);
niter = eigs2.num_iterations();
nops = eigs2.num_operations();
INFO("nconv = " << nconv);
INFO("niter = " << niter);
INFO("nops = " << nops);
REQUIRE(eigs2.info() == CompInfo::Successful);
evals = eigs2.eigenvalues();
evecs = eigs2.eigenvectors();
resid = M.selfadjointView<Eigen::Lower>() * evecs - evecs * evals.asDiagonal();
err = resid.array().abs().maxCoeff();
INFO("||AU - UD||_inf = " << err);
REQUIRE(err == Approx(0.0).margin(1e-9));
INFO("Estimated =\n " << evals);
diff = (true_evals.head(k) - evals).array().abs().maxCoeff();
INFO("diff = " << diff);
REQUIRE(diff == Approx(0.0).margin(1e-9));
}
TEST_CASE("Example #144, (n, k, m) = (20, 3, 6)", "[example_144]")
{
std::srand(123);
int n = 20;
int k = 3;
int m = 6;
run_test(n, k, m);
}
TEST_CASE("Example #144, (n, k, m) = (20, 5, 12)", "[example_144]")
{
std::srand(123);
int n = 20;
int k = 5;
int m = 12;
run_test(n, k, m);
}
TEST_CASE("Example #144, (n, k, m) = (20, 6, 12)", "[example_144]")
{
std::srand(123);
int n = 20;
int k = 6;
int m = 12;
run_test(n, k, m);
}
|
