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/GenEigsSolver.h>
8 #include <Spectra/MatOp/DenseGenMatProd.h>
9 #include <Spectra/MatOp/SparseGenMatProd.h>
10
11 using namespace Spectra;
12
13 #include "catch.hpp"
14
15 using Matrix = Eigen::MatrixXd;
16 using Vector = Eigen::VectorXd;
17 using ComplexMatrix = Eigen::MatrixXcd;
18 using ComplexVector = Eigen::VectorXcd;
19 using SpMatrix = Eigen::SparseMatrix<double>;
20
21 // Generate random sparse matrix
gen_sparse_data(int n,double prob=0.5)22 SpMatrix gen_sparse_data(int n, double prob = 0.5)
23 {
24 SpMatrix mat(n, n);
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 < n; i++)
29 {
30 for (int j = 0; j < n; j++)
31 {
32 if (distr(gen) < prob)
33 mat.insert(i, j) = distr(gen) - 0.5;
34 }
35 }
36 return mat;
37 }
38
39 template <typename MatType, typename Solver>
run_test(const MatType & mat,Solver & eigs,SortRule selection,bool allow_fail=false)40 void run_test(const MatType& mat, Solver& eigs, SortRule selection, bool allow_fail = false)
41 {
42 eigs.init();
43 // maxit = 300 to reduce running time for failed cases
44 int nconv = eigs.compute(selection, 300);
45 int niter = eigs.num_iterations();
46 int nops = eigs.num_operations();
47
48 if (allow_fail && eigs.info() != CompInfo::Successful)
49 {
50 WARN("FAILED on this test");
51 std::cout << "nconv = " << nconv << std::endl;
52 std::cout << "niter = " << niter << std::endl;
53 std::cout << "nops = " << nops << std::endl;
54 return;
55 }
56 else
57 {
58 INFO("nconv = " << nconv);
59 INFO("niter = " << niter);
60 INFO("nops = " << nops);
61 REQUIRE(eigs.info() == CompInfo::Successful);
62 }
63
64 ComplexVector evals = eigs.eigenvalues();
65 ComplexMatrix evecs = eigs.eigenvectors();
66
67 ComplexMatrix resid = mat * evecs - evecs * evals.asDiagonal();
68 const double err = resid.array().abs().maxCoeff();
69
70 INFO("||AU - UD||_inf = " << err);
71 REQUIRE(err == Approx(0.0).margin(1e-9));
72 }
73
74 template <typename MatType>
run_test_sets(const MatType & A,int k,int m)75 void run_test_sets(const MatType& A, int k, int m)
76 {
77 constexpr bool is_dense = std::is_same<MatType, Matrix>::value;
78 using DenseOp = DenseGenMatProd<double>;
79 using SparseOp = SparseGenMatProd<double>;
80 using OpType = typename std::conditional<is_dense, DenseOp, SparseOp>::type;
81
82 OpType op(A);
83 GenEigsSolver<OpType> eigs(op, k, m);
84
85 SECTION("Largest Magnitude")
86 {
87 run_test(A, eigs, SortRule::LargestMagn);
88 }
89 SECTION("Largest Real Part")
90 {
91 run_test(A, eigs, SortRule::LargestReal);
92 }
93 SECTION("Largest Imaginary Part")
94 {
95 run_test(A, eigs, SortRule::LargestImag);
96 }
97 SECTION("Smallest Magnitude")
98 {
99 run_test(A, eigs, SortRule::SmallestMagn, true);
100 }
101 SECTION("Smallest Real Part")
102 {
103 run_test(A, eigs, SortRule::SmallestReal);
104 }
105 SECTION("Smallest Imaginary Part")
106 {
107 run_test(A, eigs, SortRule::SmallestImag, true);
108 }
109 }
110
111 TEST_CASE("Eigensolver of general real matrix [10x10]", "[eigs_gen]")
112 {
113 std::srand(123);
114
115 const Matrix A = Matrix::Random(10, 10);
116 int k = 3;
117 int m = 6;
118
119 run_test_sets(A, k, m);
120 }
121
122 TEST_CASE("Eigensolver of general real matrix [100x100]", "[eigs_gen]")
123 {
124 std::srand(123);
125
126 const Matrix A = Matrix::Random(100, 100);
127 int k = 10;
128 int m = 30;
129
130 run_test_sets(A, k, m);
131 }
132
133 TEST_CASE("Eigensolver of general real matrix [1000x1000]", "[eigs_gen]")
134 {
135 std::srand(123);
136
137 const Matrix A = Matrix::Random(1000, 1000);
138 int k = 20;
139 int m = 50;
140
141 run_test_sets(A, k, m);
142 }
143
144 TEST_CASE("Eigensolver of sparse real matrix [10x10]", "[eigs_gen]")
145 {
146 std::srand(123);
147
148 const SpMatrix A = gen_sparse_data(10, 0.5);
149 int k = 3;
150 int m = 6;
151
152 run_test_sets(A, k, m);
153 }
154
155 TEST_CASE("Eigensolver of sparse real matrix [100x100]", "[eigs_gen]")
156 {
157 std::srand(123);
158
159 const SpMatrix A = gen_sparse_data(100, 0.1);
160 int k = 10;
161 int m = 30;
162
163 run_test_sets(A, k, m);
164 }
165
166 TEST_CASE("Eigensolver of sparse real matrix [1000x1000]", "[eigs_gen]")
167 {
168 std::srand(123);
169
170 const SpMatrix A = gen_sparse_data(1000, 0.01);
171 int k = 20;
172 int m = 50;
173
174 run_test_sets(A, k, m);
175 }
176