| 67 | |
| 68 | template <typename matrix_type> |
| 69 | void test_lu ( const matrix_type& m) |
| 70 | { |
| 71 | typedef typename matrix_type::type type; |
| 72 | const type eps = 10*max(abs(m))*sqrt(std::numeric_limits<type>::epsilon()); |
| 73 | dlog << LDEBUG << "test_lu(): " << m.nr() << " x " << m.nc() << " eps: " << eps; |
| 74 | print_spinner(); |
| 75 | |
| 76 | |
| 77 | lu_decomposition<matrix_type> test(m); |
| 78 | |
| 79 | DLIB_TEST(test.is_square() == (m.nr() == m.nc())); |
| 80 | |
| 81 | DLIB_TEST(test.nr() == m.nr()); |
| 82 | DLIB_TEST(test.nc() == m.nc()); |
| 83 | |
| 84 | dlog << LDEBUG << "m.nr(): " << m.nr() << " m.nc(): " << m.nc(); |
| 85 | |
| 86 | type temp; |
| 87 | DLIB_TEST_MSG( (temp= max(abs(test.get_l()*test.get_u() - rowm(m,test.get_pivot())))) < eps,temp); |
| 88 | |
| 89 | if (test.is_square()) |
| 90 | { |
| 91 | // none of the matrices we should be passing in to test_lu() should be singular. |
| 92 | DLIB_TEST_MSG (abs(test.det()) > eps/100, "det: " << test.det() ); |
| 93 | dlog << LDEBUG << "big det: " << test.det(); |
| 94 | |
| 95 | DLIB_TEST(test.is_singular() == false); |
| 96 | |
| 97 | matrix<type> m2; |
| 98 | matrix<type,0,1> col; |
| 99 | |
| 100 | m2 = identity_matrix<type>(m.nr()); |
| 101 | DLIB_TEST_MSG(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2))); |
| 102 | m2 = randmat<type>(m.nr(),5); |
| 103 | DLIB_TEST_MSG(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2))); |
| 104 | m2 = randmat<type>(m.nr(),1); |
| 105 | DLIB_TEST_MSG(equal(m*test.solve(m2), m2,eps),max(abs(m*test.solve(m2)- m2))); |
| 106 | col = randmat<type>(m.nr(),1); |
| 107 | DLIB_TEST_MSG(equal(m*test.solve(col), col,eps),max(abs(m*test.solve(m2)- m2))); |
| 108 | |
| 109 | // now make us a singular matrix |
| 110 | if (m.nr() > 1) |
| 111 | { |
| 112 | matrix<type> sm(m); |
| 113 | set_colm(sm,0) = colm(sm,1); |
| 114 | |
| 115 | lu_decomposition<matrix_type> test2(sm); |
| 116 | DLIB_TEST_MSG( (temp= max(abs(test2.get_l()*test2.get_u() - rowm(sm,test2.get_pivot())))) < eps,temp); |
| 117 | |
| 118 | // these checks are only accurate for small matrices |
| 119 | if (test2.nr() < 100) |
| 120 | { |
| 121 | DLIB_TEST_MSG(test2.is_singular() == true,"det: " << test2.det()); |
| 122 | DLIB_TEST_MSG(abs(test2.det()) < eps,"det: " << test2.det()); |
| 123 | } |
| 124 | |
| 125 | } |
| 126 | } |