eigensolver_selfadjoint.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #include "main.h"
00027 #include <limits>
00028 #include <Eigen/Eigenvalues>
00029 
00030 #ifdef HAS_GSL
00031 #include "gsl_helper.h"
00032 #endif
00033 
00034 template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
00035 {
00036   typedef typename MatrixType::Index Index;
00037   /* this test covers the following files:
00038      EigenSolver.h, SelfAdjointEigenSolver.h (and indirectly: Tridiagonalization.h)
00039   */
00040   Index rows = m.rows();
00041   Index cols = m.cols();
00042 
00043   typedef typename MatrixType::Scalar Scalar;
00044   typedef typename NumTraits<Scalar>::Real RealScalar;
00045   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00046   typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
00047   typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
00048 
00049   RealScalar largerEps = 10*test_precision<RealScalar>();
00050 
00051   MatrixType a = MatrixType::Random(rows,cols);
00052   MatrixType a1 = MatrixType::Random(rows,cols);
00053   MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
00054   symmA.template triangularView<StrictlyUpper>().setZero();
00055 
00056   MatrixType b = MatrixType::Random(rows,cols);
00057   MatrixType b1 = MatrixType::Random(rows,cols);
00058   MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1;
00059   symmB.template triangularView<StrictlyUpper>().setZero();
00060 
00061   SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);
00062   // generalized eigen pb
00063   GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB);
00064 
00065   #ifdef HAS_GSL
00066   if (internal::is_same<RealScalar,double>::value)
00067   {
00068     // restore symmA and symmB.
00069     symmA = MatrixType(symmA.template selfadjointView<Lower>());
00070     symmB = MatrixType(symmB.template selfadjointView<Lower>());
00071     typedef GslTraits<Scalar> Gsl;
00072     typename Gsl::Matrix gEvec=0, gSymmA=0, gSymmB=0;
00073     typename GslTraits<RealScalar>::Vector gEval=0;
00074     RealVectorType _eval;
00075     MatrixType _evec;
00076     convert<MatrixType>(symmA, gSymmA);
00077     convert<MatrixType>(symmB, gSymmB);
00078     convert<MatrixType>(symmA, gEvec);
00079     gEval = GslTraits<RealScalar>::createVector(rows);
00080 
00081     Gsl::eigen_symm(gSymmA, gEval, gEvec);
00082     convert(gEval, _eval);
00083     convert(gEvec, _evec);
00084 
00085     // test gsl itself !
00086     VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal(), largerEps));
00087 
00088     // compare with eigen
00089     VERIFY_IS_APPROX(_eval, eiSymm.eigenvalues());
00090     VERIFY_IS_APPROX(_evec.cwiseAbs(), eiSymm.eigenvectors().cwiseAbs());
00091 
00092     // generalized pb
00093     Gsl::eigen_symm_gen(gSymmA, gSymmB, gEval, gEvec);
00094     convert(gEval, _eval);
00095     convert(gEvec, _evec);
00096     // test GSL itself:
00097     VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps));
00098 
00099     // compare with eigen
00100     MatrixType normalized_eivec = eiSymmGen.eigenvectors()*eiSymmGen.eigenvectors().colwise().norm().asDiagonal().inverse();
00101     VERIFY_IS_APPROX(_eval, eiSymmGen.eigenvalues());
00102     VERIFY_IS_APPROX(_evec.cwiseAbs(), normalized_eivec.cwiseAbs());
00103 
00104     Gsl::free(gSymmA);
00105     Gsl::free(gSymmB);
00106     GslTraits<RealScalar>::free(gEval);
00107     Gsl::free(gEvec);
00108   }
00109   #endif
00110 
00111   VERIFY_IS_EQUAL(eiSymm.info(), Success);
00112   VERIFY((symmA.template selfadjointView<Lower>() * eiSymm.eigenvectors()).isApprox(
00113           eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
00114   VERIFY_IS_APPROX(symmA.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
00115 
00116   SelfAdjointEigenSolver<MatrixType> eiSymmNoEivecs(symmA, false);
00117   VERIFY_IS_EQUAL(eiSymmNoEivecs.info(), Success);
00118   VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues());
00119 
00120   // generalized eigen problem Ax = lBx
00121   eiSymmGen.compute(symmA, symmB,Ax_lBx);
00122   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
00123   VERIFY((symmA.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox(
00124           symmB.template selfadjointView<Lower>() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
00125 
00126   // generalized eigen problem BAx = lx
00127   eiSymmGen.compute(symmA, symmB,BAx_lx);
00128   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
00129   VERIFY((symmB.template selfadjointView<Lower>() * (symmA.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
00130          (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
00131 
00132   // generalized eigen problem ABx = lx
00133   eiSymmGen.compute(symmA, symmB,ABx_lx);
00134   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
00135   VERIFY((symmA.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
00136          (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
00137 
00138 
00139   MatrixType sqrtSymmA = eiSymm.operatorSqrt();
00140   VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA);
00141   VERIFY_IS_APPROX(sqrtSymmA, symmA.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt());
00142 
00143   MatrixType id = MatrixType::Identity(rows, cols);
00144   VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
00145 
00146   SelfAdjointEigenSolver<MatrixType> eiSymmUninitialized;
00147   VERIFY_RAISES_ASSERT(eiSymmUninitialized.info());
00148   VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvalues());
00149   VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
00150   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
00151   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
00152 
00153   eiSymmUninitialized.compute(symmA, false);
00154   VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
00155   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
00156   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
00157 
00158   // test Tridiagonalization's methods
00159   Tridiagonalization<MatrixType> tridiag(symmA);
00160   // FIXME tridiag.matrixQ().adjoint() does not work
00161   VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
00162   
00163   if (rows > 1)
00164   {
00165     // Test matrix with NaN
00166     symmA(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
00167     SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmA);
00168     VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence);
00169   }
00170 }
00171 
00172 void test_eigensolver_selfadjoint()
00173 {
00174   for(int i = 0; i < g_repeat; i++) {
00175     // very important to test a 3x3 matrix since we provide a special path for it
00176     CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) );
00177     CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
00178     CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(10,10)) );
00179     CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(19,19)) );
00180     CALL_SUBTEST_5( selfadjointeigensolver(MatrixXcd(17,17)) );
00181     
00182     CALL_SUBTEST_9( selfadjointeigensolver(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(17,17)) );
00183 
00184     // some trivial but implementation-wise tricky cases
00185     CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(1,1)) );
00186     CALL_SUBTEST_4( selfadjointeigensolver(MatrixXd(2,2)) );
00187     CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) );
00188     CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) );
00189   }
00190 
00191   // Test problem size constructors
00192   CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf>(10));
00193   CALL_SUBTEST_8(Tridiagonalization<MatrixXf>(10));
00194 }
00195 


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