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> void lu_non_invertible() 15 { 16 typedef typename MatrixType::Index Index; 17 typedef typename MatrixType::Scalar Scalar; 18 typedef typename MatrixType::RealScalar RealScalar; 19 /* this test covers the following files: 20 LU.h 21 */ 22 Index rows, cols, cols2; 23 if(MatrixType::RowsAtCompileTime==Dynamic) 24 { 25 rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); 26 } 27 else 28 { 29 rows = MatrixType::RowsAtCompileTime; 30 } 31 if(MatrixType::ColsAtCompileTime==Dynamic) 32 { 33 cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); 34 cols2 = internal::random<int>(2,EIGEN_TEST_MAX_SIZE); 35 } 36 else 37 { 38 cols2 = cols = MatrixType::ColsAtCompileTime; 39 } 40 41 enum { 42 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 43 ColsAtCompileTime = MatrixType::ColsAtCompileTime 44 }; 45 typedef typename internal::kernel_retval_base<FullPivLU<MatrixType> >::ReturnType KernelMatrixType; 46 typedef typename internal::image_retval_base<FullPivLU<MatrixType> >::ReturnType ImageMatrixType; 47 typedef Matrix<typename MatrixType::Scalar, ColsAtCompileTime, ColsAtCompileTime> 48 CMatrixType; 49 typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, RowsAtCompileTime> 50 RMatrixType; 51 52 Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1); 53 54 // The image of the zero matrix should consist of a single (zero) column vector 55 VERIFY((MatrixType::Zero(rows,cols).fullPivLu().image(MatrixType::Zero(rows,cols)).cols() == 1)); 56 57 MatrixType m1(rows, cols), m3(rows, cols2); 58 CMatrixType m2(cols, cols2); 59 createRandomPIMatrixOfRank(rank, rows, cols, m1); 60 61 FullPivLU<MatrixType> lu; 62 63 // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank 64 // of singular values are either 0 or 1. 65 // So it's not clear at all that the epsilon should play any role there. 66 lu.setThreshold(RealScalar(0.01)); 67 lu.compute(m1); 68 69 MatrixType u(rows,cols); 70 u = lu.matrixLU().template triangularView<Upper>(); 71 RMatrixType l = RMatrixType::Identity(rows,rows); 72 l.block(0,0,rows,(std::min)(rows,cols)).template triangularView<StrictlyLower>() 73 = lu.matrixLU().block(0,0,rows,(std::min)(rows,cols)); 74 75 VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u); 76 77 KernelMatrixType m1kernel = lu.kernel(); 78 ImageMatrixType m1image = lu.image(m1); 79 80 VERIFY_IS_APPROX(m1, lu.reconstructedMatrix()); 81 VERIFY(rank == lu.rank()); 82 VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); 83 VERIFY(!lu.isInjective()); 84 VERIFY(!lu.isInvertible()); 85 VERIFY(!lu.isSurjective()); 86 VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); 87 VERIFY(m1image.fullPivLu().rank() == rank); 88 VERIFY_IS_APPROX(m1 * m1.adjoint() * m1image, m1image); 89 90 m2 = CMatrixType::Random(cols,cols2); 91 m3 = m1*m2; 92 m2 = CMatrixType::Random(cols,cols2); 93 // test that the code, which does resize(), may be applied to an xpr 94 m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3); 95 VERIFY_IS_APPROX(m3, m1*m2); 96 } 97 98 template<typename MatrixType> void lu_invertible() 99 { 100 /* this test covers the following files: 101 LU.h 102 */ 103 typedef typename MatrixType::Scalar Scalar; 104 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 105 int size = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 106 107 MatrixType m1(size, size), m2(size, size), m3(size, size); 108 FullPivLU<MatrixType> lu; 109 lu.setThreshold(RealScalar(0.01)); 110 do { 111 m1 = MatrixType::Random(size,size); 112 lu.compute(m1); 113 } while(!lu.isInvertible()); 114 115 VERIFY_IS_APPROX(m1, lu.reconstructedMatrix()); 116 VERIFY(0 == lu.dimensionOfKernel()); 117 VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector 118 VERIFY(size == lu.rank()); 119 VERIFY(lu.isInjective()); 120 VERIFY(lu.isSurjective()); 121 VERIFY(lu.isInvertible()); 122 VERIFY(lu.image(m1).fullPivLu().isInvertible()); 123 m3 = MatrixType::Random(size,size); 124 m2 = lu.solve(m3); 125 VERIFY_IS_APPROX(m3, m1*m2); 126 VERIFY_IS_APPROX(m2, lu.inverse()*m3); 127 } 128 129 template<typename MatrixType> void lu_partial_piv() 130 { 131 /* this test covers the following files: 132 PartialPivLU.h 133 */ 134 typedef typename MatrixType::Index Index; 135 typedef typename MatrixType::Scalar Scalar; 136 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 137 Index rows = internal::random<Index>(1,4); 138 Index cols = rows; 139 140 MatrixType m1(cols, rows); 141 m1.setRandom(); 142 PartialPivLU<MatrixType> plu(m1); 143 144 VERIFY_IS_APPROX(m1, plu.reconstructedMatrix()); 145 } 146 147 template<typename MatrixType> void lu_verify_assert() 148 { 149 MatrixType tmp; 150 151 FullPivLU<MatrixType> lu; 152 VERIFY_RAISES_ASSERT(lu.matrixLU()) 153 VERIFY_RAISES_ASSERT(lu.permutationP()) 154 VERIFY_RAISES_ASSERT(lu.permutationQ()) 155 VERIFY_RAISES_ASSERT(lu.kernel()) 156 VERIFY_RAISES_ASSERT(lu.image(tmp)) 157 VERIFY_RAISES_ASSERT(lu.solve(tmp)) 158 VERIFY_RAISES_ASSERT(lu.determinant()) 159 VERIFY_RAISES_ASSERT(lu.rank()) 160 VERIFY_RAISES_ASSERT(lu.dimensionOfKernel()) 161 VERIFY_RAISES_ASSERT(lu.isInjective()) 162 VERIFY_RAISES_ASSERT(lu.isSurjective()) 163 VERIFY_RAISES_ASSERT(lu.isInvertible()) 164 VERIFY_RAISES_ASSERT(lu.inverse()) 165 166 PartialPivLU<MatrixType> plu; 167 VERIFY_RAISES_ASSERT(plu.matrixLU()) 168 VERIFY_RAISES_ASSERT(plu.permutationP()) 169 VERIFY_RAISES_ASSERT(plu.solve(tmp)) 170 VERIFY_RAISES_ASSERT(plu.determinant()) 171 VERIFY_RAISES_ASSERT(plu.inverse()) 172 } 173 174 void test_lu() 175 { 176 for(int i = 0; i < g_repeat; i++) { 177 CALL_SUBTEST_1( lu_non_invertible<Matrix3f>() ); 178 CALL_SUBTEST_1( lu_verify_assert<Matrix3f>() ); 179 180 CALL_SUBTEST_2( (lu_non_invertible<Matrix<double, 4, 6> >()) ); 181 CALL_SUBTEST_2( (lu_verify_assert<Matrix<double, 4, 6> >()) ); 182 183 CALL_SUBTEST_3( lu_non_invertible<MatrixXf>() ); 184 CALL_SUBTEST_3( lu_invertible<MatrixXf>() ); 185 CALL_SUBTEST_3( lu_verify_assert<MatrixXf>() ); 186 187 CALL_SUBTEST_4( lu_non_invertible<MatrixXd>() ); 188 CALL_SUBTEST_4( lu_invertible<MatrixXd>() ); 189 CALL_SUBTEST_4( lu_partial_piv<MatrixXd>() ); 190 CALL_SUBTEST_4( lu_verify_assert<MatrixXd>() ); 191 192 CALL_SUBTEST_5( lu_non_invertible<MatrixXcf>() ); 193 CALL_SUBTEST_5( lu_invertible<MatrixXcf>() ); 194 CALL_SUBTEST_5( lu_verify_assert<MatrixXcf>() ); 195 196 CALL_SUBTEST_6( lu_non_invertible<MatrixXcd>() ); 197 CALL_SUBTEST_6( lu_invertible<MatrixXcd>() ); 198 CALL_SUBTEST_6( lu_partial_piv<MatrixXcd>() ); 199 CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() ); 200 201 CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() )); 202 203 // Test problem size constructors 204 CALL_SUBTEST_9( PartialPivLU<MatrixXf>(10) ); 205 CALL_SUBTEST_9( FullPivLU<MatrixXf>(10, 20); ); 206 } 207 } 208
