| 62 | |
| 63 | template <typename matrix_type, typename U> |
| 64 | void test_eigenvalue_impl ( const matrix_type& m, const eigenvalue_decomposition<U>& test ) |
| 65 | { |
| 66 | typedef typename matrix_type::type type; |
| 67 | const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon()); |
| 68 | dlog << LDEBUG << "test_eigenvalue(): " << m.nr() << " x " << m.nc() << " eps: " << eps; |
| 69 | print_spinner(); |
| 70 | |
| 71 | |
| 72 | DLIB_TEST(test.dim() == m.nr()); |
| 73 | |
| 74 | // make sure all the various ways of asking for the eigenvalues are actually returning a |
| 75 | // consistent set of eigenvalues. |
| 76 | DLIB_TEST(equal(real(test.get_eigenvalues()), test.get_real_eigenvalues(), eps)); |
| 77 | DLIB_TEST(equal(imag(test.get_eigenvalues()), test.get_imag_eigenvalues(), eps)); |
| 78 | DLIB_TEST(equal(real(diag(test.get_d())), test.get_real_eigenvalues(), eps)); |
| 79 | DLIB_TEST(equal(imag(diag(test.get_d())), test.get_imag_eigenvalues(), eps)); |
| 80 | |
| 81 | matrix<type> eig1 ( real_eigenvalues(m)); |
| 82 | matrix<type> eig2 ( test.get_real_eigenvalues()); |
| 83 | sort(&eig1(0), &eig1(0) + eig1.size()); |
| 84 | sort(&eig2(0), &eig2(0) + eig2.size()); |
| 85 | DLIB_TEST(max(abs(eig1 - eig2)) < eps); |
| 86 | |
| 87 | const matrix<type> V = test.get_pseudo_v(); |
| 88 | const matrix<type> D = test.get_pseudo_d(); |
| 89 | const matrix<complex<type> > CV = test.get_v(); |
| 90 | const matrix<complex<type> > CD = test.get_d(); |
| 91 | const matrix<complex<type> > CM = complex_matrix(m, uniform_matrix<type>(m.nr(),m.nc(),0)); |
| 92 | |
| 93 | DLIB_TEST(V.nr() == test.dim()); |
| 94 | DLIB_TEST(V.nc() == test.dim()); |
| 95 | DLIB_TEST(D.nr() == test.dim()); |
| 96 | DLIB_TEST(D.nc() == test.dim()); |
| 97 | |
| 98 | // CD is a diagonal matrix |
| 99 | DLIB_TEST(diagm(diag(CD)) == CD); |
| 100 | |
| 101 | // verify that these things are actually eigenvalues and eigenvectors of m |
| 102 | DLIB_TEST_MSG(max(abs(m*V - V*D)) < eps, max(abs(m*V - V*D)) << " " << eps); |
| 103 | DLIB_TEST(max(norm(CM*CV - CV*CD)) < eps); |
| 104 | |
| 105 | // if m is a symmetric matrix |
| 106 | if (max(abs(m-trans(m))) < 1e-5) |
| 107 | { |
| 108 | dlog << LTRACE << "m is symmetric"; |
| 109 | // there aren't any imaginary eigenvalues |
| 110 | DLIB_TEST(max(abs(test.get_imag_eigenvalues())) < eps); |
| 111 | DLIB_TEST(diagm(diag(D)) == D); |
| 112 | |
| 113 | // only check the determinant against the eigenvalues for small matrices |
| 114 | // because for huge ones the determinant might be so big it overflows a floating point number. |
| 115 | if (m.nr() < 50) |
| 116 | { |
| 117 | const type mdet = det(m); |
| 118 | DLIB_TEST_MSG(std::abs(prod(test.get_real_eigenvalues()) - mdet) < std::abs(mdet)*sqrt(std::numeric_limits<type>::epsilon()), |
| 119 | std::abs(prod(test.get_real_eigenvalues()) - mdet) <<" eps: " << std::abs(mdet)*sqrt(std::numeric_limits<type>::epsilon()) |
| 120 | << " mdet: "<< mdet << " prod(eig): " << prod(test.get_real_eigenvalues()) |
| 121 | ); |
no test coverage detected