gtsam
3rdparty
Eigen
test
eigensolver_generalized_real.cpp
Go to the documentation of this file.
1
// This file is part of Eigen, a lightweight C++ template library
2
// for linear algebra.
3
//
4
// Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
5
//
6
// This Source Code Form is subject to the terms of the Mozilla
7
// Public License v. 2.0. If a copy of the MPL was not distributed
8
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10
#define EIGEN_RUNTIME_NO_MALLOC
11
#include "
main.h
"
12
#include <limits>
13
#include <Eigen/Eigenvalues>
14
#include <Eigen/LU>
15
16
template
<
typename
MatrixType>
void
generalized_eigensolver_real
(
const
MatrixType
&
m
)
17
{
18
/* this test covers the following files:
19
GeneralizedEigenSolver.h
20
*/
21
Index
rows
=
m
.rows();
22
Index
cols
=
m
.cols();
23
24
typedef
typename
MatrixType::Scalar
Scalar
;
25
typedef
std::complex<Scalar> ComplexScalar;
26
typedef
Matrix<Scalar, MatrixType::RowsAtCompileTime, 1>
VectorType
;
27
28
MatrixType
a
= MatrixType::Random(
rows
,
cols
);
29
MatrixType
b
= MatrixType::Random(
rows
,
cols
);
30
MatrixType
a1
= MatrixType::Random(
rows
,
cols
);
31
MatrixType
b1
= MatrixType::Random(
rows
,
cols
);
32
MatrixType
spdA =
a
.adjoint() *
a
+
a1
.adjoint() *
a1
;
33
MatrixType
spdB =
b
.adjoint() *
b
+
b1
.adjoint() *
b1
;
34
35
// lets compare to GeneralizedSelfAdjointEigenSolver
36
{
37
GeneralizedSelfAdjointEigenSolver<MatrixType>
symmEig(spdA, spdB);
38
GeneralizedEigenSolver<MatrixType>
eig
(spdA, spdB);
39
40
VERIFY_IS_EQUAL
(
eig
.eigenvalues().imag().cwiseAbs().maxCoeff(), 0);
41
42
VectorType
realEigenvalues =
eig
.eigenvalues().real();
43
std::sort(realEigenvalues.data(), realEigenvalues.data()+realEigenvalues.size());
44
VERIFY_IS_APPROX
(realEigenvalues, symmEig.
eigenvalues
());
45
46
// check eigenvectors
47
typename
GeneralizedEigenSolver<MatrixType>::EigenvectorsType
D
=
eig
.eigenvalues().asDiagonal();
48
typename
GeneralizedEigenSolver<MatrixType>::EigenvectorsType
V
=
eig
.eigenvectors();
49
VERIFY_IS_APPROX
(spdA*
V
, spdB*
V
*
D
);
50
}
51
52
// non symmetric case:
53
{
54
GeneralizedEigenSolver<MatrixType>
eig
(
rows
);
55
// TODO enable full-prealocation of required memory, this probably requires an in-place mode for HessenbergDecomposition
56
//Eigen::internal::set_is_malloc_allowed(false);
57
eig
.compute(
a
,
b
);
58
//Eigen::internal::set_is_malloc_allowed(true);
59
for
(
Index
k=0; k<
cols
; ++k)
60
{
61
Matrix<ComplexScalar,Dynamic,Dynamic>
tmp = (
eig
.betas()(k)*
a
).template cast<ComplexScalar>() -
eig
.alphas()(k)*
b
;
62
if
(tmp.size()>1 && tmp.norm()>(
std::numeric_limits<Scalar>::min
)())
63
tmp /= tmp.norm();
64
VERIFY_IS_MUCH_SMALLER_THAN
(
std::abs
(tmp.determinant()),
Scalar
(1) );
65
}
66
// check eigenvectors
67
typename
GeneralizedEigenSolver<MatrixType>::EigenvectorsType
D
=
eig
.eigenvalues().asDiagonal();
68
typename
GeneralizedEigenSolver<MatrixType>::EigenvectorsType
V
=
eig
.eigenvectors();
69
VERIFY_IS_APPROX
(
a
*
V
,
b
*
V
*
D
);
70
}
71
72
// regression test for bug 1098
73
{
74
GeneralizedSelfAdjointEigenSolver<MatrixType>
eig1(
a
.adjoint() *
a
,
b
.adjoint() *
b
);
75
eig1.
compute
(
a
.adjoint() *
a
,
b
.adjoint() *
b
);
76
GeneralizedEigenSolver<MatrixType>
eig2(
a
.adjoint() *
a
,
b
.adjoint() *
b
);
77
eig2.
compute
(
a
.adjoint() *
a
,
b
.adjoint() *
b
);
78
}
79
80
// check without eigenvectors
81
{
82
GeneralizedEigenSolver<MatrixType>
eig1(spdA, spdB,
true
);
83
GeneralizedEigenSolver<MatrixType>
eig2(spdA, spdB,
false
);
84
VERIFY_IS_APPROX
(eig1.
eigenvalues
(), eig2.
eigenvalues
());
85
}
86
}
87
88
EIGEN_DECLARE_TEST
(eigensolver_generalized_real)
89
{
90
for
(
int
i
= 0;
i
<
g_repeat
;
i
++) {
91
int
s
= 0;
92
CALL_SUBTEST_1
(
generalized_eigensolver_real
(Matrix4f()) );
93
s
= internal::random<int>(1,
EIGEN_TEST_MAX_SIZE
/4);
94
CALL_SUBTEST_2
(
generalized_eigensolver_real
(MatrixXd(
s
,
s
)) );
95
96
// some trivial but implementation-wise special cases
97
CALL_SUBTEST_2
(
generalized_eigensolver_real
(MatrixXd(1,1)) );
98
CALL_SUBTEST_2
(
generalized_eigensolver_real
(MatrixXd(2,2)) );
99
CALL_SUBTEST_3
(
generalized_eigensolver_real
(
Matrix<double,1,1>
()) );
100
CALL_SUBTEST_4
(
generalized_eigensolver_real
(Matrix2d()) );
101
TEST_SET_BUT_UNUSED_VARIABLE
(
s
)
102
}
103
}
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition:
main.h:390
simple_graph::b1
Vector2 b1(2, -1)
D
MatrixXcd D
Definition:
EigenSolver_EigenSolver_MatrixType.cpp:14
s
RealScalar s
Definition:
level1_cplx_impl.h:126
Eigen::GeneralizedEigenSolver
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
Definition:
GeneralizedEigenSolver.h:58
MatrixType
MatrixXf MatrixType
Definition:
benchmark-blocking-sizes.cpp:52
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition:
main.h:386
b
Scalar * b
Definition:
benchVecAdd.cpp:17
Eigen::GeneralizedEigenSolver::compute
GeneralizedEigenSolver & compute(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Computes generalized eigendecomposition of given matrix.
Definition:
GeneralizedEigenSolver.h:287
Eigen::SelfAdjointEigenSolver::eigenvalues
const EIGEN_DEVICE_FUNC RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition:
SelfAdjointEigenSolver.h:300
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(eigensolver_generalized_real)
Definition:
eigensolver_generalized_real.cpp:88
rows
int rows
Definition:
Tutorial_commainit_02.cpp:1
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition:
split_test_helper.h:22
align_3::a1
Point2 a1
Definition:
testPose2.cpp:769
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition:
split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition:
split_test_helper.h:4
generalized_eigensolver_real
void generalized_eigensolver_real(const MatrixType &m)
Definition:
eigensolver_generalized_real.cpp:16
Eigen::g_repeat
static int g_repeat
Definition:
main.h:169
m
Matrix3f m
Definition:
AngleAxis_mimic_euler.cpp:1
eig
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition:
split_test_helper.h:10
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition:
integer_types.cpp:15
a
ArrayXXi a
Definition:
Array_initializer_list_23_cxx11.cpp:1
main.h
Eigen::GeneralizedSelfAdjointEigenSolver
Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem.
Definition:
GeneralizedSelfAdjointEigenSolver.h:48
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition:
boostmultiprec.cpp:16
min
#define min(a, b)
Definition:
datatypes.h:19
V
MatrixXcd V
Definition:
EigenSolver_EigenSolver_MatrixType.cpp:15
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition:
3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
abs
#define abs(x)
Definition:
datatypes.h:17
VectorType
Definition:
FFTW.cpp:65
Eigen::GeneralizedEigenSolver::eigenvalues
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Definition:
GeneralizedEigenSolver.h:202
Eigen::GeneralizedSelfAdjointEigenSolver::compute
GeneralizedSelfAdjointEigenSolver & compute(const MatrixType &matA, const MatrixType &matB, int options=ComputeEigenvectors|Ax_lBx)
Computes generalized eigendecomposition of given matrix pencil.
Definition:
GeneralizedSelfAdjointEigenSolver.h:163
cols
int cols
Definition:
Tutorial_commainit_02.cpp:1
i
int i
Definition:
BiCGSTAB_step_by_step.cpp:9
TEST_SET_BUT_UNUSED_VARIABLE
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition:
main.h:121
Scalar
SCALAR Scalar
Definition:
bench_gemm.cpp:46
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition:
Meta.h:74
gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:29