eigen2_eigensolver.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "main.h"
00026 #include <Eigen/QR>
00027 
00028 #ifdef HAS_GSL
00029 #include "gsl_helper.h"
00030 #endif
00031 
00032 template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
00033 {
00034   /* this test covers the following files:
00035      EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h)
00036   */
00037   int rows = m.rows();
00038   int cols = m.cols();
00039 
00040   typedef typename MatrixType::Scalar Scalar;
00041   typedef typename NumTraits<Scalar>::Real RealScalar;
00042   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00043   typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
00044   typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
00045 
00046   RealScalar largerEps = 10*test_precision<RealScalar>();
00047 
00048   MatrixType a = MatrixType::Random(rows,cols);
00049   MatrixType a1 = MatrixType::Random(rows,cols);
00050   MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
00051 
00052   MatrixType b = MatrixType::Random(rows,cols);
00053   MatrixType b1 = MatrixType::Random(rows,cols);
00054   MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
00055 
00056   SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);
00057   // generalized eigen pb
00058   SelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB);
00059 
00060   #ifdef HAS_GSL
00061   if (ei_is_same_type<RealScalar,double>::ret)
00062   {
00063     typedef GslTraits<Scalar> Gsl;
00064     typename Gsl::Matrix gEvec=0, gSymmA=0, gSymmB=0;
00065     typename GslTraits<RealScalar>::Vector gEval=0;
00066     RealVectorType _eval;
00067     MatrixType _evec;
00068     convert<MatrixType>(symmA, gSymmA);
00069     convert<MatrixType>(symmB, gSymmB);
00070     convert<MatrixType>(symmA, gEvec);
00071     gEval = GslTraits<RealScalar>::createVector(rows);
00072 
00073     Gsl::eigen_symm(gSymmA, gEval, gEvec);
00074     convert(gEval, _eval);
00075     convert(gEvec, _evec);
00076 
00077     // test gsl itself !
00078     VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal(), largerEps));
00079 
00080     // compare with eigen
00081     VERIFY_IS_APPROX(_eval, eiSymm.eigenvalues());
00082     VERIFY_IS_APPROX(_evec.cwise().abs(), eiSymm.eigenvectors().cwise().abs());
00083 
00084     // generalized pb
00085     Gsl::eigen_symm_gen(gSymmA, gSymmB, gEval, gEvec);
00086     convert(gEval, _eval);
00087     convert(gEvec, _evec);
00088     // test GSL itself:
00089     VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps));
00090 
00091     // compare with eigen
00092     MatrixType normalized_eivec = eiSymmGen.eigenvectors()*eiSymmGen.eigenvectors().colwise().norm().asDiagonal().inverse();
00093     VERIFY_IS_APPROX(_eval, eiSymmGen.eigenvalues());
00094     VERIFY_IS_APPROX(_evec.cwiseAbs(), normalized_eivec.cwiseAbs());
00095 
00096     Gsl::free(gSymmA);
00097     Gsl::free(gSymmB);
00098     GslTraits<RealScalar>::free(gEval);
00099     Gsl::free(gEvec);
00100   }
00101   #endif
00102 
00103   VERIFY((symmA * eiSymm.eigenvectors()).isApprox(
00104           eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
00105 
00106   // generalized eigen problem Ax = lBx
00107   VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
00108           symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
00109 
00110   MatrixType sqrtSymmA = eiSymm.operatorSqrt();
00111   VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA);
00112   VERIFY_IS_APPROX(sqrtSymmA, symmA*eiSymm.operatorInverseSqrt());
00113 }
00114 
00115 template<typename MatrixType> void eigensolver(const MatrixType& m)
00116 {
00117   /* this test covers the following files:
00118      EigenSolver.h
00119   */
00120   int rows = m.rows();
00121   int cols = m.cols();
00122 
00123   typedef typename MatrixType::Scalar Scalar;
00124   typedef typename NumTraits<Scalar>::Real RealScalar;
00125   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00126   typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
00127   typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
00128 
00129   // RealScalar largerEps = 10*test_precision<RealScalar>();
00130 
00131   MatrixType a = MatrixType::Random(rows,cols);
00132   MatrixType a1 = MatrixType::Random(rows,cols);
00133   MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
00134 
00135   EigenSolver<MatrixType> ei0(symmA);
00136   VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix());
00137   VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
00138     (ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
00139 
00140   EigenSolver<MatrixType> ei1(a);
00141   VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
00142   VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
00143                    ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
00144 
00145 }
00146 
00147 void test_eigen2_eigensolver()
00148 {
00149   for(int i = 0; i < g_repeat; i++) {
00150     // very important to test a 3x3 matrix since we provide a special path for it
00151     CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) );
00152     CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
00153     CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(7,7)) );
00154     CALL_SUBTEST_4( selfadjointeigensolver(MatrixXcd(5,5)) );
00155     CALL_SUBTEST_5( selfadjointeigensolver(MatrixXd(19,19)) );
00156 
00157     CALL_SUBTEST_6( eigensolver(Matrix4f()) );
00158     CALL_SUBTEST_5( eigensolver(MatrixXd(17,17)) );
00159   }
00160 }
00161 


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:38