18bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)// This file is part of Eigen, a lightweight C++ template library 28bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)// for linear algebra. Eigen itself is part of the KDE project. 38bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)// 48bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> 58bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)// 68bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)// This Source Code Form is subject to the terms of the Mozilla 78bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)// Public License v. 2.0. If a copy of the MPL was not distributed 8a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles)// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 923730a6e56a168d1879203e4b3819bb36e3d8f1fTorne (Richard Coles) 108bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)#include "main.h" 118bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)#include <Eigen/LU> 128bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 138bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)template<typename Derived> 141e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles)void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m) 151e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles){ 168bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) typedef typename Derived::RealScalar RealScalar; 17f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) for(int a = 0; a < 3*(m.rows()+m.cols()); a++) 188bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) { 198bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) RealScalar d = Eigen::ei_random<RealScalar>(-1,1); 201e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number 211e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) int j; 221e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) do { 238bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) j = Eigen::ei_random<int>(0,m.rows()-1); 248bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) } while (i==j); // j is another one (must be different) 258bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) m.row(i) += d * m.row(j); 268bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 278bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number 281e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) do { 298bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) j = Eigen::ei_random<int>(0,m.cols()-1); 308bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) } while (i==j); // j is another one (must be different) 318bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) m.col(i) += d * m.col(j); 328bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) } 338bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)} 348bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 358bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)template<typename MatrixType> void lu_non_invertible() 368bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles){ 378bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) /* this test covers the following files: 388bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) LU.h 398bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) */ 408bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) // NOTE there seems to be a problem with too small sizes -- could easily lie in the doSomeRankPreservingOperations function 418bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200); 428bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) int rank = ei_random<int>(1, std::min(rows, cols)-1); 438bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 448bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1); 458bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) m1 = MatrixType::Random(rows,cols); 468bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) if(rows <= cols) 478bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) for(int i = rank; i < rows; i++) m1.row(i).setZero(); 488bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) else 498bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) for(int i = rank; i < cols; i++) m1.col(i).setZero(); 508bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) doSomeRankPreservingOperations(m1); 518bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 528bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) LU<MatrixType> lu(m1); 538bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) typename LU<MatrixType>::KernelResultType m1kernel = lu.kernel(); 548bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) typename LU<MatrixType>::ImageResultType m1image = lu.image(); 558bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 56a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) VERIFY(rank == lu.rank()); 578bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); 588bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) VERIFY(!lu.isInjective()); 598bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) VERIFY(!lu.isInvertible()); 608bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) VERIFY(lu.isSurjective() == (lu.rank() == rows)); 618bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); 628bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) VERIFY(m1image.lu().rank() == rank); 638bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); 641e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) sidebyside << m1, m1image; 651e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) VERIFY(sidebyside.lu().rank() == rank); 661e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) m2 = MatrixType::Random(cols,cols2); 671e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) m3 = m1*m2; 681e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) m2 = MatrixType::Random(cols,cols2); 691e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) lu.solve(m3, &m2); 701e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) VERIFY_IS_APPROX(m3, m1*m2); 71a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) /* solve now always returns true 72a1401311d1ab56c4ed0a474bd38c108f75cb0cd9Torne (Richard Coles) m3 = MatrixType::Random(rows,cols2); 731e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) VERIFY(!lu.solve(m3, &m2)); 741e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles) */ 751e9bf3e0803691d0a228da41fc608347b6db4340Torne (Richard Coles)} 768bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 778bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles)template<typename MatrixType> void lu_invertible() 788bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles){ 798bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) /* this test covers the following files: 808bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) LU.h 818bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) */ 828bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 838bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) int size = ei_random<int>(10,200); 848bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 858bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) MatrixType m1(size, size), m2(size, size), m3(size, size); 868bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) m1 = MatrixType::Random(size,size); 878bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) 888bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) if (ei_is_same_type<RealScalar,float>::ret) 89f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) { 90f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) // let's build a matrix more stable to inverse 91f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) MatrixType a = MatrixType::Random(size,size*2); 92f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) m1 += a * a.adjoint(); 93f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) } 94f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) 95f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) LU<MatrixType> lu(m1); 96f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) VERIFY(0 == lu.dimensionOfKernel()); 97f2477e01787aa58f445919b809d89e252beef54fTorne (Richard Coles) VERIFY(size == lu.rank()); 988bcbed890bc3ce4d7a057a8f32cab53fa534672eTorne (Richard Coles) VERIFY(lu.isInjective()); 99 VERIFY(lu.isSurjective()); 100 VERIFY(lu.isInvertible()); 101 VERIFY(lu.image().lu().isInvertible()); 102 m3 = MatrixType::Random(size,size); 103 lu.solve(m3, &m2); 104 VERIFY_IS_APPROX(m3, m1*m2); 105 VERIFY_IS_APPROX(m2, lu.inverse()*m3); 106 m3 = MatrixType::Random(size,size); 107 VERIFY(lu.solve(m3, &m2)); 108} 109 110void test_eigen2_lu() 111{ 112 for(int i = 0; i < g_repeat; i++) { 113 CALL_SUBTEST_1( lu_non_invertible<MatrixXf>() ); 114 CALL_SUBTEST_2( lu_non_invertible<MatrixXd>() ); 115 CALL_SUBTEST_3( lu_non_invertible<MatrixXcf>() ); 116 CALL_SUBTEST_4( lu_non_invertible<MatrixXcd>() ); 117 CALL_SUBTEST_1( lu_invertible<MatrixXf>() ); 118 CALL_SUBTEST_2( lu_invertible<MatrixXd>() ); 119 CALL_SUBTEST_3( lu_invertible<MatrixXcf>() ); 120 CALL_SUBTEST_4( lu_invertible<MatrixXcd>() ); 121 } 122} 123