1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2009 Benoit Jacob <jacob.benoit.1 (at) gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #include "main.h" 11 #include <Eigen/LU> 12 using namespace std; 13 14 template<typename MatrixType> 15 typename MatrixType::RealScalar matrix_l1_norm(const MatrixType& m) { 16 return m.cwiseAbs().colwise().sum().maxCoeff(); 17 } 18 19 template<typename MatrixType> void lu_non_invertible() 20 { 21 typedef typename MatrixType::Index Index; 22 typedef typename MatrixType::RealScalar RealScalar; 23 /* this test covers the following files: 24 LU.h 25 */ 26 Index rows, cols, cols2; 27 if(MatrixType::RowsAtCompileTime==Dynamic) 28 { 29 rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); 30 } 31 else 32 { 33 rows = MatrixType::RowsAtCompileTime; 34 } 35 if(MatrixType::ColsAtCompileTime==Dynamic) 36 { 37 cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); 38 cols2 = internal::random<int>(2,EIGEN_TEST_MAX_SIZE); 39 } 40 else 41 { 42 cols2 = cols = MatrixType::ColsAtCompileTime; 43 } 44 45 enum { 46 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 47 ColsAtCompileTime = MatrixType::ColsAtCompileTime 48 }; 49 typedef typename internal::kernel_retval_base<FullPivLU<MatrixType> >::ReturnType KernelMatrixType; 50 typedef typename internal::image_retval_base<FullPivLU<MatrixType> >::ReturnType ImageMatrixType; 51 typedef Matrix<typename MatrixType::Scalar, ColsAtCompileTime, ColsAtCompileTime> 52 CMatrixType; 53 typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, RowsAtCompileTime> 54 RMatrixType; 55 56 Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1); 57 58 // The image of the zero matrix should consist of a single (zero) column vector 59 VERIFY((MatrixType::Zero(rows,cols).fullPivLu().image(MatrixType::Zero(rows,cols)).cols() == 1)); 60 61 MatrixType m1(rows, cols), m3(rows, cols2); 62 CMatrixType m2(cols, cols2); 63 createRandomPIMatrixOfRank(rank, rows, cols, m1); 64 65 FullPivLU<MatrixType> lu; 66 67 // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank 68 // of singular values are either 0 or 1. 69 // So it's not clear at all that the epsilon should play any role there. 70 lu.setThreshold(RealScalar(0.01)); 71 lu.compute(m1); 72 73 MatrixType u(rows,cols); 74 u = lu.matrixLU().template triangularView<Upper>(); 75 RMatrixType l = RMatrixType::Identity(rows,rows); 76 l.block(0,0,rows,(std::min)(rows,cols)).template triangularView<StrictlyLower>() 77 = lu.matrixLU().block(0,0,rows,(std::min)(rows,cols)); 78 79 VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u); 80 81 KernelMatrixType m1kernel = lu.kernel(); 82 ImageMatrixType m1image = lu.image(m1); 83 84 VERIFY_IS_APPROX(m1, lu.reconstructedMatrix()); 85 VERIFY(rank == lu.rank()); 86 VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); 87 VERIFY(!lu.isInjective()); 88 VERIFY(!lu.isInvertible()); 89 VERIFY(!lu.isSurjective()); 90 VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); 91 VERIFY(m1image.fullPivLu().rank() == rank); 92 VERIFY_IS_APPROX(m1 * m1.adjoint() * m1image, m1image); 93 94 m2 = CMatrixType::Random(cols,cols2); 95 m3 = m1*m2; 96 m2 = CMatrixType::Random(cols,cols2); 97 // test that the code, which does resize(), may be applied to an xpr 98 m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3); 99 VERIFY_IS_APPROX(m3, m1*m2); 100 101 // test solve with transposed 102 m3 = MatrixType::Random(rows,cols2); 103 m2 = m1.transpose()*m3; 104 m3 = MatrixType::Random(rows,cols2); 105 lu.template _solve_impl_transposed<false>(m2, m3); 106 VERIFY_IS_APPROX(m2, m1.transpose()*m3); 107 m3 = MatrixType::Random(rows,cols2); 108 m3 = lu.transpose().solve(m2); 109 VERIFY_IS_APPROX(m2, m1.transpose()*m3); 110 111 // test solve with conjugate transposed 112 m3 = MatrixType::Random(rows,cols2); 113 m2 = m1.adjoint()*m3; 114 m3 = MatrixType::Random(rows,cols2); 115 lu.template _solve_impl_transposed<true>(m2, m3); 116 VERIFY_IS_APPROX(m2, m1.adjoint()*m3); 117 m3 = MatrixType::Random(rows,cols2); 118 m3 = lu.adjoint().solve(m2); 119 VERIFY_IS_APPROX(m2, m1.adjoint()*m3); 120 } 121 122 template<typename MatrixType> void lu_invertible() 123 { 124 /* this test covers the following files: 125 LU.h 126 */ 127 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 128 Index size = MatrixType::RowsAtCompileTime; 129 if( size==Dynamic) 130 size = internal::random<Index>(1,EIGEN_TEST_MAX_SIZE); 131 132 MatrixType m1(size, size), m2(size, size), m3(size, size); 133 FullPivLU<MatrixType> lu; 134 lu.setThreshold(RealScalar(0.01)); 135 do { 136 m1 = MatrixType::Random(size,size); 137 lu.compute(m1); 138 } while(!lu.isInvertible()); 139 140 VERIFY_IS_APPROX(m1, lu.reconstructedMatrix()); 141 VERIFY(0 == lu.dimensionOfKernel()); 142 VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector 143 VERIFY(size == lu.rank()); 144 VERIFY(lu.isInjective()); 145 VERIFY(lu.isSurjective()); 146 VERIFY(lu.isInvertible()); 147 VERIFY(lu.image(m1).fullPivLu().isInvertible()); 148 m3 = MatrixType::Random(size,size); 149 m2 = lu.solve(m3); 150 VERIFY_IS_APPROX(m3, m1*m2); 151 MatrixType m1_inverse = lu.inverse(); 152 VERIFY_IS_APPROX(m2, m1_inverse*m3); 153 154 RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); 155 const RealScalar rcond_est = lu.rcond(); 156 // Verify that the estimated condition number is within a factor of 10 of the 157 // truth. 158 VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); 159 160 // test solve with transposed 161 lu.template _solve_impl_transposed<false>(m3, m2); 162 VERIFY_IS_APPROX(m3, m1.transpose()*m2); 163 m3 = MatrixType::Random(size,size); 164 m3 = lu.transpose().solve(m2); 165 VERIFY_IS_APPROX(m2, m1.transpose()*m3); 166 167 // test solve with conjugate transposed 168 lu.template _solve_impl_transposed<true>(m3, m2); 169 VERIFY_IS_APPROX(m3, m1.adjoint()*m2); 170 m3 = MatrixType::Random(size,size); 171 m3 = lu.adjoint().solve(m2); 172 VERIFY_IS_APPROX(m2, m1.adjoint()*m3); 173 174 // Regression test for Bug 302 175 MatrixType m4 = MatrixType::Random(size,size); 176 VERIFY_IS_APPROX(lu.solve(m3*m4), lu.solve(m3)*m4); 177 } 178 179 template<typename MatrixType> void lu_partial_piv() 180 { 181 /* this test covers the following files: 182 PartialPivLU.h 183 */ 184 typedef typename MatrixType::Index Index; 185 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 186 Index size = internal::random<Index>(1,4); 187 188 MatrixType m1(size, size), m2(size, size), m3(size, size); 189 m1.setRandom(); 190 PartialPivLU<MatrixType> plu(m1); 191 192 VERIFY_IS_APPROX(m1, plu.reconstructedMatrix()); 193 194 m3 = MatrixType::Random(size,size); 195 m2 = plu.solve(m3); 196 VERIFY_IS_APPROX(m3, m1*m2); 197 MatrixType m1_inverse = plu.inverse(); 198 VERIFY_IS_APPROX(m2, m1_inverse*m3); 199 200 RealScalar rcond = (RealScalar(1) / matrix_l1_norm(m1)) / matrix_l1_norm(m1_inverse); 201 const RealScalar rcond_est = plu.rcond(); 202 // Verify that the estimate is within a factor of 10 of the truth. 203 VERIFY(rcond_est > rcond / 10 && rcond_est < rcond * 10); 204 205 // test solve with transposed 206 plu.template _solve_impl_transposed<false>(m3, m2); 207 VERIFY_IS_APPROX(m3, m1.transpose()*m2); 208 m3 = MatrixType::Random(size,size); 209 m3 = plu.transpose().solve(m2); 210 VERIFY_IS_APPROX(m2, m1.transpose()*m3); 211 212 // test solve with conjugate transposed 213 plu.template _solve_impl_transposed<true>(m3, m2); 214 VERIFY_IS_APPROX(m3, m1.adjoint()*m2); 215 m3 = MatrixType::Random(size,size); 216 m3 = plu.adjoint().solve(m2); 217 VERIFY_IS_APPROX(m2, m1.adjoint()*m3); 218 } 219 220 template<typename MatrixType> void lu_verify_assert() 221 { 222 MatrixType tmp; 223 224 FullPivLU<MatrixType> lu; 225 VERIFY_RAISES_ASSERT(lu.matrixLU()) 226 VERIFY_RAISES_ASSERT(lu.permutationP()) 227 VERIFY_RAISES_ASSERT(lu.permutationQ()) 228 VERIFY_RAISES_ASSERT(lu.kernel()) 229 VERIFY_RAISES_ASSERT(lu.image(tmp)) 230 VERIFY_RAISES_ASSERT(lu.solve(tmp)) 231 VERIFY_RAISES_ASSERT(lu.determinant()) 232 VERIFY_RAISES_ASSERT(lu.rank()) 233 VERIFY_RAISES_ASSERT(lu.dimensionOfKernel()) 234 VERIFY_RAISES_ASSERT(lu.isInjective()) 235 VERIFY_RAISES_ASSERT(lu.isSurjective()) 236 VERIFY_RAISES_ASSERT(lu.isInvertible()) 237 VERIFY_RAISES_ASSERT(lu.inverse()) 238 239 PartialPivLU<MatrixType> plu; 240 VERIFY_RAISES_ASSERT(plu.matrixLU()) 241 VERIFY_RAISES_ASSERT(plu.permutationP()) 242 VERIFY_RAISES_ASSERT(plu.solve(tmp)) 243 VERIFY_RAISES_ASSERT(plu.determinant()) 244 VERIFY_RAISES_ASSERT(plu.inverse()) 245 } 246 247 void test_lu() 248 { 249 for(int i = 0; i < g_repeat; i++) { 250 CALL_SUBTEST_1( lu_non_invertible<Matrix3f>() ); 251 CALL_SUBTEST_1( lu_invertible<Matrix3f>() ); 252 CALL_SUBTEST_1( lu_verify_assert<Matrix3f>() ); 253 254 CALL_SUBTEST_2( (lu_non_invertible<Matrix<double, 4, 6> >()) ); 255 CALL_SUBTEST_2( (lu_verify_assert<Matrix<double, 4, 6> >()) ); 256 257 CALL_SUBTEST_3( lu_non_invertible<MatrixXf>() ); 258 CALL_SUBTEST_3( lu_invertible<MatrixXf>() ); 259 CALL_SUBTEST_3( lu_verify_assert<MatrixXf>() ); 260 261 CALL_SUBTEST_4( lu_non_invertible<MatrixXd>() ); 262 CALL_SUBTEST_4( lu_invertible<MatrixXd>() ); 263 CALL_SUBTEST_4( lu_partial_piv<MatrixXd>() ); 264 CALL_SUBTEST_4( lu_verify_assert<MatrixXd>() ); 265 266 CALL_SUBTEST_5( lu_non_invertible<MatrixXcf>() ); 267 CALL_SUBTEST_5( lu_invertible<MatrixXcf>() ); 268 CALL_SUBTEST_5( lu_verify_assert<MatrixXcf>() ); 269 270 CALL_SUBTEST_6( lu_non_invertible<MatrixXcd>() ); 271 CALL_SUBTEST_6( lu_invertible<MatrixXcd>() ); 272 CALL_SUBTEST_6( lu_partial_piv<MatrixXcd>() ); 273 CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() ); 274 275 CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() )); 276 277 // Test problem size constructors 278 CALL_SUBTEST_9( PartialPivLU<MatrixXf>(10) ); 279 CALL_SUBTEST_9( FullPivLU<MatrixXf>(10, 20); ); 280 } 281 } 282
