Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "main.h"
00027 #include <limits>
00028 #include <Eigen/Eigenvalues>
00029 #include <Eigen/LU>
00030
00031
00032
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
00052
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
00074
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
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
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
00124 CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(10));
00125 }