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 (at) 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 17 template<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 100 template<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 132 void 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
