1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra. Eigen itself is part of the KDE project.
3//
4// Copyright (C) 2008 Gael Guennebaud <g.gael@free.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#include "main.h"
11#include <Eigen/QR>
12
13#ifdef HAS_GSL
14#include "gsl_helper.h"
15#endif
16
17template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
18{
19  /* this test covers the following files:
20     EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h)
21  */
22  int rows = m.rows();
23  int cols = m.cols();
24
25  typedef typename MatrixType::Scalar Scalar;
26  typedef typename NumTraits<Scalar>::Real RealScalar;
27  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
28  typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
29  typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
30
31  RealScalar largerEps = 10*test_precision<RealScalar>();
32
33  MatrixType a = MatrixType::Random(rows,cols);
34  MatrixType a1 = MatrixType::Random(rows,cols);
35  MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
36
37  MatrixType b = MatrixType::Random(rows,cols);
38  MatrixType b1 = MatrixType::Random(rows,cols);
39  MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
40
41  SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);
42  // generalized eigen pb
43  SelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB);
44
45  #ifdef HAS_GSL
46  if (ei_is_same_type<RealScalar,double>::ret)
47  {
48    typedef GslTraits<Scalar> Gsl;
49    typename Gsl::Matrix gEvec=0, gSymmA=0, gSymmB=0;
50    typename GslTraits<RealScalar>::Vector gEval=0;
51    RealVectorType _eval;
52    MatrixType _evec;
53    convert<MatrixType>(symmA, gSymmA);
54    convert<MatrixType>(symmB, gSymmB);
55    convert<MatrixType>(symmA, gEvec);
56    gEval = GslTraits<RealScalar>::createVector(rows);
57
58    Gsl::eigen_symm(gSymmA, gEval, gEvec);
59    convert(gEval, _eval);
60    convert(gEvec, _evec);
61
62    // test gsl itself !
63    VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal(), largerEps));
64
65    // compare with eigen
66    VERIFY_IS_APPROX(_eval, eiSymm.eigenvalues());
67    VERIFY_IS_APPROX(_evec.cwise().abs(), eiSymm.eigenvectors().cwise().abs());
68
69    // generalized pb
70    Gsl::eigen_symm_gen(gSymmA, gSymmB, gEval, gEvec);
71    convert(gEval, _eval);
72    convert(gEvec, _evec);
73    // test GSL itself:
74    VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps));
75
76    // compare with eigen
77    MatrixType normalized_eivec = eiSymmGen.eigenvectors()*eiSymmGen.eigenvectors().colwise().norm().asDiagonal().inverse();
78    VERIFY_IS_APPROX(_eval, eiSymmGen.eigenvalues());
79    VERIFY_IS_APPROX(_evec.cwiseAbs(), normalized_eivec.cwiseAbs());
80
81    Gsl::free(gSymmA);
82    Gsl::free(gSymmB);
83    GslTraits<RealScalar>::free(gEval);
84    Gsl::free(gEvec);
85  }
86  #endif
87
88  VERIFY((symmA * eiSymm.eigenvectors()).isApprox(
89          eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
90
91  // generalized eigen problem Ax = lBx
92  VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
93          symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
94
95  MatrixType sqrtSymmA = eiSymm.operatorSqrt();
96  VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA);
97  VERIFY_IS_APPROX(sqrtSymmA, symmA*eiSymm.operatorInverseSqrt());
98}
99
100template<typename MatrixType> void eigensolver(const MatrixType& m)
101{
102  /* this test covers the following files:
103     EigenSolver.h
104  */
105  int rows = m.rows();
106  int cols = m.cols();
107
108  typedef typename MatrixType::Scalar Scalar;
109  typedef typename NumTraits<Scalar>::Real RealScalar;
110  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
111  typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
112  typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
113
114  // RealScalar largerEps = 10*test_precision<RealScalar>();
115
116  MatrixType a = MatrixType::Random(rows,cols);
117  MatrixType a1 = MatrixType::Random(rows,cols);
118  MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
119
120  EigenSolver<MatrixType> ei0(symmA);
121  VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix());
122  VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
123    (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
124
125  EigenSolver<MatrixType> ei1(a);
126  VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
127  VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
128                   ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
129
130}
131
132void test_eigen2_eigensolver()
133{
134  for(int i = 0; i < g_repeat; i++) {
135    // very important to test a 3x3 matrix since we provide a special path for it
136    CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) );
137    CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
138    CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(7,7)) );
139    CALL_SUBTEST_4( selfadjointeigensolver(MatrixXcd(5,5)) );
140    CALL_SUBTEST_5( selfadjointeigensolver(MatrixXd(19,19)) );
141
142    CALL_SUBTEST_6( eigensolver(Matrix4f()) );
143    CALL_SUBTEST_5( eigensolver(MatrixXd(17,17)) );
144  }
145}
146
147