1 #include <Eigen/Core>
2 #include <Eigen/SparseCore>
3 #include <iostream>
4 #include <type_traits>
5 #include <random> // Requires C++ 11
6
7 #include <Spectra/SymGEigsSolver.h>
8 #include <Spectra/MatOp/DenseSymMatProd.h>
9 #include <Spectra/MatOp/DenseCholesky.h>
10 #include <Spectra/MatOp/SparseSymMatProd.h>
11 #include <Spectra/MatOp/SparseCholesky.h>
12
13 using namespace Spectra;
14
15 #include "catch.hpp"
16
17 using Matrix = Eigen::MatrixXd;
18 using Vector = Eigen::VectorXd;
19 using SpMatrix = Eigen::SparseMatrix<double>;
20
21 // Generate random sparse matrix
sprand(int size,double prob=0.5)22 SpMatrix sprand(int size, double prob = 0.5)
23 {
24 SpMatrix mat(size, size);
25 std::default_random_engine gen;
26 gen.seed(0);
27 std::uniform_real_distribution<double> distr(0.0, 1.0);
28 for (int i = 0; i < size; i++)
29 {
30 for (int j = 0; j < size; j++)
31 {
32 if (distr(gen) < prob)
33 mat.insert(i, j) = distr(gen) - 0.5;
34 }
35 }
36 return mat;
37 }
38
39 // Generate data for testing
gen_dense_data(int n,Matrix & A,Matrix & B)40 void gen_dense_data(int n, Matrix& A, Matrix& B)
41 {
42 Matrix M = Eigen::MatrixXd::Random(n, n);
43 A = M + M.transpose();
44 B = M.transpose() * M;
45 // To make sure B is positive definite
46 B.diagonal() += Eigen::VectorXd::Random(n).cwiseAbs();
47 }
48
gen_sparse_data(int n,SpMatrix & A,SpMatrix & B,double prob=0.1)49 void gen_sparse_data(int n, SpMatrix& A, SpMatrix& B, double prob = 0.1)
50 {
51 // Eigen solver only uses the lower triangle of A,
52 // so we don't need to make A symmetric here.
53 A = sprand(n, prob);
54 B = A.transpose() * A;
55 // To make sure B is positive definite
56 for (int i = 0; i < n; i++)
57 B.coeffRef(i, i) += 0.1;
58 }
59
60 template <typename MatType, typename Solver>
run_test(const MatType & A,const MatType & B,Solver & eigs,SortRule selection,bool allow_fail=false)61 void run_test(const MatType& A, const MatType& B, Solver& eigs, SortRule selection, bool allow_fail = false)
62 {
63 eigs.init();
64 // maxit = 300 to reduce running time for failed cases
65 int nconv = eigs.compute(selection, 300);
66 int niter = eigs.num_iterations();
67 int nops = eigs.num_operations();
68
69 if (allow_fail && eigs.info() != CompInfo::Successful)
70 {
71 WARN("FAILED on this test");
72 std::cout << "nconv = " << nconv << std::endl;
73 std::cout << "niter = " << niter << std::endl;
74 std::cout << "nops = " << nops << std::endl;
75 return;
76 }
77 else
78 {
79 INFO("nconv = " << nconv);
80 INFO("niter = " << niter);
81 INFO("nops = " << nops);
82 REQUIRE(eigs.info() == CompInfo::Successful);
83 }
84
85 Vector evals = eigs.eigenvalues();
86 Matrix evecs = eigs.eigenvectors();
87
88 Matrix resid = A.template selfadjointView<Eigen::Lower>() * evecs -
89 B.template selfadjointView<Eigen::Lower>() * evecs * evals.asDiagonal();
90 const double err = resid.array().abs().maxCoeff();
91
92 INFO("||AU - BUD||_inf = " << err);
93 REQUIRE(err == Approx(0.0).margin(1e-9));
94 }
95
96 template <typename MatType>
run_test_sets(const MatType & A,const MatType & B,int k,int m)97 void run_test_sets(const MatType& A, const MatType& B, int k, int m)
98 {
99 constexpr bool is_dense = std::is_same<MatType, Matrix>::value;
100 using DenseOp = DenseSymMatProd<double>;
101 using SparseOp = SparseSymMatProd<double>;
102 using OpType = typename std::conditional<is_dense, DenseOp, SparseOp>::type;
103
104 using DenseBOp = DenseCholesky<double>;
105 using SparseBOp = SparseCholesky<double>;
106 using BOpType = typename std::conditional<is_dense, DenseBOp, SparseBOp>::type;
107
108 OpType op(A);
109 BOpType Bop(B);
110 // Make sure B is positive definite and the decomposition is successful
111 REQUIRE(Bop.info() == CompInfo::Successful);
112
113 SymGEigsSolver<OpType, BOpType, GEigsMode::Cholesky> eigs(op, Bop, k, m);
114
115 SECTION("Largest Magnitude")
116 {
117 run_test(A, B, eigs, SortRule::LargestMagn);
118 }
119 SECTION("Largest Value")
120 {
121 run_test(A, B, eigs, SortRule::LargestAlge);
122 }
123 SECTION("Smallest Magnitude")
124 {
125 run_test(A, B, eigs, SortRule::SmallestMagn, true);
126 }
127 SECTION("Smallest Value")
128 {
129 run_test(A, B, eigs, SortRule::SmallestAlge);
130 }
131 SECTION("Both Ends")
132 {
133 run_test(A, B, eigs, SortRule::BothEnds);
134 }
135 }
136
137 TEST_CASE("Generalized eigensolver of symmetric real matrix [10x10]", "[geigs_sym]")
138 {
139 std::srand(123);
140
141 Matrix A, B;
142 gen_dense_data(10, A, B);
143 int k = 3;
144 int m = 6;
145
146 run_test_sets(A, B, k, m);
147 }
148
149 TEST_CASE("Generalized eigensolver of symmetric real matrix [100x100]", "[geigs_sym]")
150 {
151 std::srand(123);
152
153 Matrix A, B;
154 gen_dense_data(100, A, B);
155 int k = 10;
156 int m = 20;
157
158 run_test_sets(A, B, k, m);
159 }
160
161 TEST_CASE("Generalized eigensolver of symmetric real matrix [1000x1000]", "[geigs_sym]")
162 {
163 std::srand(123);
164
165 Matrix A, B;
166 gen_dense_data(1000, A, B);
167 int k = 20;
168 int m = 50;
169
170 run_test_sets(A, B, k, m);
171 }
172
173 TEST_CASE("Generalized eigensolver of sparse symmetric real matrix [10x10]", "[geigs_sym]")
174 {
175 std::srand(123);
176
177 // Eigen solver only uses the lower triangle
178 SpMatrix A, B;
179 gen_sparse_data(10, A, B, 0.5);
180 int k = 3;
181 int m = 6;
182
183 run_test_sets(A, B, k, m);
184 }
185
186 TEST_CASE("Generalized eigensolver of sparse symmetric real matrix [100x100]", "[geigs_sym]")
187 {
188 std::srand(123);
189
190 // Eigen solver only uses the lower triangle
191 SpMatrix A, B;
192 gen_sparse_data(100, A, B, 0.1);
193 int k = 10;
194 int m = 20;
195
196 run_test_sets(A, B, k, m);
197 }
198
199 TEST_CASE("Generalized eigensolver of sparse symmetric real matrix [1000x1000]", "[geigs_sym]")
200 {
201 std::srand(123);
202
203 // Eigen solver only uses the lower triangle
204 SpMatrix A, B;
205 gen_sparse_data(1000, A, B, 0.01);
206 int k = 20;
207 int m = 50;
208
209 run_test_sets(A, B, k, m);
210 }
211