1a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles)// This file is part of Eigen, a lightweight C++ template library 258537e28ecd584eab876aee8be7156509866d23aTorne (Richard Coles)// for linear algebra. 358537e28ecd584eab876aee8be7156509866d23aTorne (Richard Coles)// 458537e28ecd584eab876aee8be7156509866d23aTorne (Richard Coles)// Copyright (C) 2008-2009 Benoit Jacob <jacob.benoit.1@gmail.com> 5a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles)// 6a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles)// This Source Code Form is subject to the terms of the Mozilla 758537e28ecd584eab876aee8be7156509866d23aTorne (Richard Coles)// Public License v. 2.0. If a copy of the MPL was not distributed 85f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 95f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles) 1058537e28ecd584eab876aee8be7156509866d23aTorne (Richard Coles)#include "main.h" 11a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles)#include <Eigen/LU> 125f1c94371a64b3196d4be9466099bb892df9b88eTorne (Richard Coles)using namespace std; 1358537e28ecd584eab876aee8be7156509866d23aTorne (Richard Coles) 1458537e28ecd584eab876aee8be7156509866d23aTorne (Richard Coles)template<typename MatrixType> void lu_non_invertible() 156e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles){ 166e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) typedef typename MatrixType::Index Index; 176e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) typedef typename MatrixType::RealScalar RealScalar; 186e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) /* this test covers the following files: 196e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) LU.h 206e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) */ 216e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) Index rows, cols, cols2; 226e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) if(MatrixType::RowsAtCompileTime==Dynamic) 236e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) { 246e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); 256e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) } 266e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) else 276e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) { 286e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) rows = MatrixType::RowsAtCompileTime; 296e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) } 306e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) if(MatrixType::ColsAtCompileTime==Dynamic) 316e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) { 326e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); 336e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) cols2 = internal::random<int>(2,EIGEN_TEST_MAX_SIZE); 346e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) } 356e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) else 366e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) { 376e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) cols2 = cols = MatrixType::ColsAtCompileTime; 386e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) } 396e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) 406e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) enum { 416e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) RowsAtCompileTime = MatrixType::RowsAtCompileTime, 426e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) ColsAtCompileTime = MatrixType::ColsAtCompileTime 436e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) }; 446e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) typedef typename internal::kernel_retval_base<FullPivLU<MatrixType> >::ReturnType KernelMatrixType; 456e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) typedef typename internal::image_retval_base<FullPivLU<MatrixType> >::ReturnType ImageMatrixType; 466e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) typedef Matrix<typename MatrixType::Scalar, ColsAtCompileTime, ColsAtCompileTime> 476e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) CMatrixType; 486e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, RowsAtCompileTime> 496e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) RMatrixType; 506e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) 516e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1); 526e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) 536e8cce623b6e4fe0c9e4af605d675dd9d0338c38Torne (Richard Coles) // The image of the zero matrix should consist of a single (zero) column vector 5458537e28ecd584eab876aee8be7156509866d23aTorne (Richard Coles) VERIFY((MatrixType::Zero(rows,cols).fullPivLu().image(MatrixType::Zero(rows,cols)).cols() == 1)); 55a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) 56 MatrixType m1(rows, cols), m3(rows, cols2); 57 CMatrixType m2(cols, cols2); 58 createRandomPIMatrixOfRank(rank, rows, cols, m1); 59 60 FullPivLU<MatrixType> lu; 61 62 // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank 63 // of singular values are either 0 or 1. 64 // So it's not clear at all that the epsilon should play any role there. 65 lu.setThreshold(RealScalar(0.01)); 66 lu.compute(m1); 67 68 MatrixType u(rows,cols); 69 u = lu.matrixLU().template triangularView<Upper>(); 70 RMatrixType l = RMatrixType::Identity(rows,rows); 71 l.block(0,0,rows,(std::min)(rows,cols)).template triangularView<StrictlyLower>() 72 = lu.matrixLU().block(0,0,rows,(std::min)(rows,cols)); 73 74 VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u); 75 76 KernelMatrixType m1kernel = lu.kernel(); 77 ImageMatrixType m1image = lu.image(m1); 78 79 VERIFY_IS_APPROX(m1, lu.reconstructedMatrix()); 80 VERIFY(rank == lu.rank()); 81 VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); 82 VERIFY(!lu.isInjective()); 83 VERIFY(!lu.isInvertible()); 84 VERIFY(!lu.isSurjective()); 85 VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); 86 VERIFY(m1image.fullPivLu().rank() == rank); 87 VERIFY_IS_APPROX(m1 * m1.adjoint() * m1image, m1image); 88 89 m2 = CMatrixType::Random(cols,cols2); 90 m3 = m1*m2; 91 m2 = CMatrixType::Random(cols,cols2); 92 // test that the code, which does resize(), may be applied to an xpr 93 m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3); 94 VERIFY_IS_APPROX(m3, m1*m2); 95} 96 97template<typename MatrixType> void lu_invertible() 98{ 99 /* this test covers the following files: 100 LU.h 101 */ 102 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 103 int size = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 104 105 MatrixType m1(size, size), m2(size, size), m3(size, size); 106 FullPivLU<MatrixType> lu; 107 lu.setThreshold(RealScalar(0.01)); 108 do { 109 m1 = MatrixType::Random(size,size); 110 lu.compute(m1); 111 } while(!lu.isInvertible()); 112 113 VERIFY_IS_APPROX(m1, lu.reconstructedMatrix()); 114 VERIFY(0 == lu.dimensionOfKernel()); 115 VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector 116 VERIFY(size == lu.rank()); 117 VERIFY(lu.isInjective()); 118 VERIFY(lu.isSurjective()); 119 VERIFY(lu.isInvertible()); 120 VERIFY(lu.image(m1).fullPivLu().isInvertible()); 121 m3 = MatrixType::Random(size,size); 122 m2 = lu.solve(m3); 123 VERIFY_IS_APPROX(m3, m1*m2); 124 VERIFY_IS_APPROX(m2, lu.inverse()*m3); 125} 126 127template<typename MatrixType> void lu_partial_piv() 128{ 129 /* this test covers the following files: 130 PartialPivLU.h 131 */ 132 typedef typename MatrixType::Index Index; 133 Index rows = internal::random<Index>(1,4); 134 Index cols = rows; 135 136 MatrixType m1(cols, rows); 137 m1.setRandom(); 138 PartialPivLU<MatrixType> plu(m1); 139 140 VERIFY_IS_APPROX(m1, plu.reconstructedMatrix()); 141} 142 143template<typename MatrixType> void lu_verify_assert() 144{ 145 MatrixType tmp; 146 147 FullPivLU<MatrixType> lu; 148 VERIFY_RAISES_ASSERT(lu.matrixLU()) 149 VERIFY_RAISES_ASSERT(lu.permutationP()) 150 VERIFY_RAISES_ASSERT(lu.permutationQ()) 151 VERIFY_RAISES_ASSERT(lu.kernel()) 152 VERIFY_RAISES_ASSERT(lu.image(tmp)) 153 VERIFY_RAISES_ASSERT(lu.solve(tmp)) 154 VERIFY_RAISES_ASSERT(lu.determinant()) 155 VERIFY_RAISES_ASSERT(lu.rank()) 156 VERIFY_RAISES_ASSERT(lu.dimensionOfKernel()) 157 VERIFY_RAISES_ASSERT(lu.isInjective()) 158 VERIFY_RAISES_ASSERT(lu.isSurjective()) 159 VERIFY_RAISES_ASSERT(lu.isInvertible()) 160 VERIFY_RAISES_ASSERT(lu.inverse()) 161 162 PartialPivLU<MatrixType> plu; 163 VERIFY_RAISES_ASSERT(plu.matrixLU()) 164 VERIFY_RAISES_ASSERT(plu.permutationP()) 165 VERIFY_RAISES_ASSERT(plu.solve(tmp)) 166 VERIFY_RAISES_ASSERT(plu.determinant()) 167 VERIFY_RAISES_ASSERT(plu.inverse()) 168} 169 170void test_lu() 171{ 172 for(int i = 0; i < g_repeat; i++) { 173 CALL_SUBTEST_1( lu_non_invertible<Matrix3f>() ); 174 CALL_SUBTEST_1( lu_verify_assert<Matrix3f>() ); 175 176 CALL_SUBTEST_2( (lu_non_invertible<Matrix<double, 4, 6> >()) ); 177 CALL_SUBTEST_2( (lu_verify_assert<Matrix<double, 4, 6> >()) ); 178 179 CALL_SUBTEST_3( lu_non_invertible<MatrixXf>() ); 180 CALL_SUBTEST_3( lu_invertible<MatrixXf>() ); 181 CALL_SUBTEST_3( lu_verify_assert<MatrixXf>() ); 182 183 CALL_SUBTEST_4( lu_non_invertible<MatrixXd>() ); 184 CALL_SUBTEST_4( lu_invertible<MatrixXd>() ); 185 CALL_SUBTEST_4( lu_partial_piv<MatrixXd>() ); 186 CALL_SUBTEST_4( lu_verify_assert<MatrixXd>() ); 187 188 CALL_SUBTEST_5( lu_non_invertible<MatrixXcf>() ); 189 CALL_SUBTEST_5( lu_invertible<MatrixXcf>() ); 190 CALL_SUBTEST_5( lu_verify_assert<MatrixXcf>() ); 191 192 CALL_SUBTEST_6( lu_non_invertible<MatrixXcd>() ); 193 CALL_SUBTEST_6( lu_invertible<MatrixXcd>() ); 194 CALL_SUBTEST_6( lu_partial_piv<MatrixXcd>() ); 195 CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() ); 196 197 CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() )); 198 199 // Test problem size constructors 200 CALL_SUBTEST_9( PartialPivLU<MatrixXf>(10) ); 201 CALL_SUBTEST_9( FullPivLU<MatrixXf>(10, 20); ); 202 } 203} 204