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