eigensolver_complex.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-2009 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 #include <Eigen/LU>
00030 
00031 /* Check that two column vectors are approximately equal upto permutations,
00032    by checking that the k-th power sums are equal for k = 1, ..., vec1.rows() */
00033 template<typename VectorType>
00034 void verify_is_approx_upto_permutation(const VectorType& vec1, const VectorType& vec2)
00035 {
00036   typedef typename NumTraits<typename VectorType::Scalar>::Real RealScalar;
00037 
00038   VERIFY(vec1.cols() == 1);
00039   VERIFY(vec2.cols() == 1);
00040   VERIFY(vec1.rows() == vec2.rows());
00041   for (int k = 1; k <= vec1.rows(); ++k)
00042   {
00043     VERIFY_IS_APPROX(vec1.array().pow(RealScalar(k)).sum(), vec2.array().pow(RealScalar(k)).sum());
00044   }
00045 }
00046 
00047 
00048 template<typename MatrixType> void eigensolver(const MatrixType& m)
00049 {
00050   typedef typename MatrixType::Index Index;
00051   /* this test covers the following files:
00052      ComplexEigenSolver.h, and indirectly ComplexSchur.h
00053   */
00054   Index rows = m.rows();
00055   Index cols = m.cols();
00056 
00057   typedef typename MatrixType::Scalar Scalar;
00058   typedef typename NumTraits<Scalar>::Real RealScalar;
00059   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00060   typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
00061   typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
00062 
00063   MatrixType a = MatrixType::Random(rows,cols);
00064   MatrixType symmA =  a.adjoint() * a;
00065 
00066   ComplexEigenSolver<MatrixType> ei0(symmA);
00067   VERIFY_IS_EQUAL(ei0.info(), Success);
00068   VERIFY_IS_APPROX(symmA * ei0.eigenvectors(), ei0.eigenvectors() * ei0.eigenvalues().asDiagonal());
00069 
00070   ComplexEigenSolver<MatrixType> ei1(a);
00071   VERIFY_IS_EQUAL(ei1.info(), Success);
00072   VERIFY_IS_APPROX(a * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
00073   // Note: If MatrixType is real then a.eigenvalues() uses EigenSolver and thus
00074   // another algorithm so results may differ slightly
00075   verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues());
00076 
00077   ComplexEigenSolver<MatrixType> eiNoEivecs(a, false);
00078   VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
00079   VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
00080 
00081   // Regression test for issue #66
00082   MatrixType z = MatrixType::Zero(rows,cols);
00083   ComplexEigenSolver<MatrixType> eiz(z);
00084   VERIFY((eiz.eigenvalues().cwiseEqual(0)).all());
00085 
00086   MatrixType id = MatrixType::Identity(rows, cols);
00087   VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
00088 
00089   if (rows > 1)
00090   {
00091     // Test matrix with NaN
00092     a(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
00093     ComplexEigenSolver<MatrixType> eiNaN(a);
00094     VERIFY_IS_EQUAL(eiNaN.info(), NoConvergence);
00095   }
00096 }
00097 
00098 template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
00099 {
00100   ComplexEigenSolver<MatrixType> eig;
00101   VERIFY_RAISES_ASSERT(eig.eigenvectors());
00102   VERIFY_RAISES_ASSERT(eig.eigenvalues());
00103 
00104   MatrixType a = MatrixType::Random(m.rows(),m.cols());
00105   eig.compute(a, false);
00106   VERIFY_RAISES_ASSERT(eig.eigenvectors());
00107 }
00108 
00109 void test_eigensolver_complex()
00110 {
00111   for(int i = 0; i < g_repeat; i++) {
00112     CALL_SUBTEST_1( eigensolver(Matrix4cf()) );
00113     CALL_SUBTEST_2( eigensolver(MatrixXcd(14,14)) );
00114     CALL_SUBTEST_3( eigensolver(Matrix<std::complex<float>, 1, 1>()) );
00115     CALL_SUBTEST_4( eigensolver(Matrix3f()) );
00116   }
00117 
00118   CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4cf()) );
00119   CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXcd(14,14)) );
00120   CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<std::complex<float>, 1, 1>()) );
00121   CALL_SUBTEST_4( eigensolver_verify_assert(Matrix3f()) );
00122 
00123   // Test problem size constructors
00124   CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(10));
00125 }


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