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
|
#include <iostream>
#include <random>
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Spectra/SymGEigsSolver.h>
#include <Spectra/MatOp/SparseSymMatProd.h>
#include <Spectra/MatOp/SparseRegularInverse.h>
#include "catch.hpp"
using namespace Spectra;
using Matrix = Eigen::MatrixXd;
using Vector = Eigen::VectorXd;
using SpMatrix = Eigen::SparseMatrix<double>;
// Generate random sparse matrix
SpMatrix sprand(int size, double prob = 0.5)
{
SpMatrix mat(size, size);
std::default_random_engine gen;
gen.seed(0);
std::uniform_real_distribution<double> distr(0.0, 1.0);
for (int i = 0; i < size; i++)
{
for (int j = 0; j < size; j++)
{
if (distr(gen) < prob)
mat.insert(i, j) = distr(gen) - 0.5;
}
}
return mat;
}
void gen_sparse_data(int n, SpMatrix& A, SpMatrix& B, double prob = 0.1)
{
// Eigen solver only uses the lower triangle of A,
// so we don't need to make A symmetric here.
A = sprand(n, prob);
B = A.transpose() * A;
// To make sure B is positive definite
for (int i = 0; i < n; i++)
B.coeffRef(i, i) += 0.1;
}
template <typename Solver>
void run_test(const SpMatrix& A, const SpMatrix& B, Solver& eigs, SortRule selection, bool allow_fail = false)
{
eigs.init();
// maxit = 100 to reduce running time for failed cases
int nconv = eigs.compute(selection, 100);
int niter = eigs.num_iterations();
int nops = eigs.num_operations();
if (allow_fail && eigs.info() != CompInfo::Successful)
{
WARN("FAILED on this test");
std::cout << "nconv = " << nconv << std::endl;
std::cout << "niter = " << niter << std::endl;
std::cout << "nops = " << nops << std::endl;
return;
}
else
{
INFO("nconv = " << nconv);
INFO("niter = " << niter);
INFO("nops = " << nops);
REQUIRE(eigs.info() == CompInfo::Successful);
}
Vector evals = eigs.eigenvalues();
Matrix evecs = eigs.eigenvectors();
Matrix resid = A.template selfadjointView<Eigen::Lower>() * evecs -
B.template selfadjointView<Eigen::Lower>() * evecs * evals.asDiagonal();
const double err = resid.array().abs().maxCoeff();
INFO("||AU - BUD||_inf = " << err);
REQUIRE(err == Approx(0.0).margin(1e-9));
}
void run_test_sets(const SpMatrix& A, const SpMatrix& B, int k, int m)
{
using OpType = SparseSymMatProd<double>;
using BOpType = SparseRegularInverse<double>;
OpType op(A);
BOpType Bop(B);
SymGEigsSolver<OpType, BOpType, GEigsMode::RegularInverse> eigs(op, Bop, k, m);
SECTION("Largest Magnitude")
{
run_test(A, B, eigs, SortRule::LargestMagn);
}
SECTION("Largest Value")
{
run_test(A, B, eigs, SortRule::LargestAlge);
}
SECTION("Smallest Magnitude")
{
run_test(A, B, eigs, SortRule::SmallestMagn, true);
}
SECTION("Smallest Value")
{
run_test(A, B, eigs, SortRule::SmallestAlge);
}
SECTION("Both Ends")
{
run_test(A, B, eigs, SortRule::BothEnds);
}
}
TEST_CASE("Generalized eigensolver of sparse symmetric real matrix [10x10]", "[geigs_sym]")
{
std::srand(123);
// Eigen solver only uses the lower triangle
SpMatrix A, B;
gen_sparse_data(10, A, B, 0.5);
int k = 3;
int m = 6;
run_test_sets(A, B, k, m);
}
TEST_CASE("Generalized eigensolver of sparse symmetric real matrix [100x100]", "[geigs_sym]")
{
std::srand(123);
// Eigen solver only uses the lower triangle
SpMatrix A, B;
gen_sparse_data(100, A, B, 0.1);
int k = 10;
int m = 20;
run_test_sets(A, B, k, m);
}
TEST_CASE("Generalized eigensolver of sparse symmetric real matrix [1000x1000]", "[geigs_sym]")
{
std::srand(123);
// Eigen solver only uses the lower triangle
SpMatrix A, B;
gen_sparse_data(1000, A, B, 0.01);
int k = 20;
int m = 50;
run_test_sets(A, B, k, m);
}
|
