1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. Eigen itself is part of the KDE project. 3 // 4 // Copyright (C) 2008 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 13 template<typename Derived> 14 void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m) 15 { 16 typedef typename Derived::RealScalar RealScalar; 17 for(int a = 0; a < 3*(m.rows()+m.cols()); a++) 18 { 19 RealScalar d = Eigen::ei_random<RealScalar>(-1,1); 20 int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number 21 int j; 22 do { 23 j = Eigen::ei_random<int>(0,m.rows()-1); 24 } while (i==j); // j is another one (must be different) 25 m.row(i) += d * m.row(j); 26 27 i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number 28 do { 29 j = Eigen::ei_random<int>(0,m.cols()-1); 30 } while (i==j); // j is another one (must be different) 31 m.col(i) += d * m.col(j); 32 } 33 } 34 35 template<typename MatrixType> void lu_non_invertible() 36 { 37 /* this test covers the following files: 38 LU.h 39 */ 40 // NOTE there seems to be a problem with too small sizes -- could easily lie in the doSomeRankPreservingOperations function 41 int rows = ei_random<int>(20,200), cols = ei_random<int>(20,200), cols2 = ei_random<int>(20,200); 42 int rank = ei_random<int>(1, std::min(rows, cols)-1); 43 44 MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1); 45 m1 = MatrixType::Random(rows,cols); 46 if(rows <= cols) 47 for(int i = rank; i < rows; i++) m1.row(i).setZero(); 48 else 49 for(int i = rank; i < cols; i++) m1.col(i).setZero(); 50 doSomeRankPreservingOperations(m1); 51 52 LU<MatrixType> lu(m1); 53 typename LU<MatrixType>::KernelResultType m1kernel = lu.kernel(); 54 typename LU<MatrixType>::ImageResultType m1image = lu.image(); 55 56 VERIFY(rank == lu.rank()); 57 VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); 58 VERIFY(!lu.isInjective()); 59 VERIFY(!lu.isInvertible()); 60 VERIFY(lu.isSurjective() == (lu.rank() == rows)); 61 VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); 62 VERIFY(m1image.lu().rank() == rank); 63 MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); 64 sidebyside << m1, m1image; 65 VERIFY(sidebyside.lu().rank() == rank); 66 m2 = MatrixType::Random(cols,cols2); 67 m3 = m1*m2; 68 m2 = MatrixType::Random(cols,cols2); 69 lu.solve(m3, &m2); 70 VERIFY_IS_APPROX(m3, m1*m2); 71 /* solve now always returns true 72 m3 = MatrixType::Random(rows,cols2); 73 VERIFY(!lu.solve(m3, &m2)); 74 */ 75 } 76 77 template<typename MatrixType> void lu_invertible() 78 { 79 /* this test covers the following files: 80 LU.h 81 */ 82 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 83 int size = ei_random<int>(10,200); 84 85 MatrixType m1(size, size), m2(size, size), m3(size, size); 86 m1 = MatrixType::Random(size,size); 87 88 if (ei_is_same_type<RealScalar,float>::ret) 89 { 90 // let's build a matrix more stable to inverse 91 MatrixType a = MatrixType::Random(size,size*2); 92 m1 += a * a.adjoint(); 93 } 94 95 LU<MatrixType> lu(m1); 96 VERIFY(0 == lu.dimensionOfKernel()); 97 VERIFY(size == lu.rank()); 98 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 110 void 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
