qr_fullpivoting.cpp (5600B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #include "main.h" 12 #include <Eigen/QR> 13 #include "solverbase.h" 14 15 template<typename MatrixType> void qr() 16 { 17 STATIC_CHECK(( internal::is_same<typename FullPivHouseholderQR<MatrixType>::StorageIndex,int>::value )); 18 19 static const int Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime; 20 Index max_size = EIGEN_TEST_MAX_SIZE; 21 Index min_size = numext::maxi(1,EIGEN_TEST_MAX_SIZE/10); 22 Index rows = Rows == Dynamic ? internal::random<Index>(min_size,max_size) : Rows, 23 cols = Cols == Dynamic ? internal::random<Index>(min_size,max_size) : Cols, 24 cols2 = Cols == Dynamic ? internal::random<Index>(min_size,max_size) : Cols, 25 rank = internal::random<Index>(1, (std::min)(rows, cols)-1); 26 27 typedef typename MatrixType::Scalar Scalar; 28 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; 29 MatrixType m1; 30 createRandomPIMatrixOfRank(rank,rows,cols,m1); 31 FullPivHouseholderQR<MatrixType> qr(m1); 32 VERIFY_IS_EQUAL(rank, qr.rank()); 33 VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel()); 34 VERIFY(!qr.isInjective()); 35 VERIFY(!qr.isInvertible()); 36 VERIFY(!qr.isSurjective()); 37 38 MatrixType r = qr.matrixQR(); 39 40 MatrixQType q = qr.matrixQ(); 41 VERIFY_IS_UNITARY(q); 42 43 // FIXME need better way to construct trapezoid 44 for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); 45 46 MatrixType c = qr.matrixQ() * r * qr.colsPermutation().inverse(); 47 48 VERIFY_IS_APPROX(m1, c); 49 50 // stress the ReturnByValue mechanism 51 MatrixType tmp; 52 VERIFY_IS_APPROX(tmp.noalias() = qr.matrixQ() * r, (qr.matrixQ() * r).eval()); 53 54 check_solverbase<MatrixType, MatrixType>(m1, qr, rows, cols, cols2); 55 56 { 57 MatrixType m2, m3; 58 Index size = rows; 59 do { 60 m1 = MatrixType::Random(size,size); 61 qr.compute(m1); 62 } while(!qr.isInvertible()); 63 MatrixType m1_inv = qr.inverse(); 64 m3 = m1 * MatrixType::Random(size,cols2); 65 m2 = qr.solve(m3); 66 VERIFY_IS_APPROX(m2, m1_inv*m3); 67 } 68 } 69 70 template<typename MatrixType> void qr_invertible() 71 { 72 using std::log; 73 using std::abs; 74 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 75 typedef typename MatrixType::Scalar Scalar; 76 77 Index max_size = numext::mini(50,EIGEN_TEST_MAX_SIZE); 78 Index min_size = numext::maxi(1,EIGEN_TEST_MAX_SIZE/10); 79 Index size = internal::random<Index>(min_size,max_size); 80 81 MatrixType m1(size, size), m2(size, size), m3(size, size); 82 m1 = MatrixType::Random(size,size); 83 84 if (internal::is_same<RealScalar,float>::value) 85 { 86 // let's build a matrix more stable to inverse 87 MatrixType a = MatrixType::Random(size,size*2); 88 m1 += a * a.adjoint(); 89 } 90 91 FullPivHouseholderQR<MatrixType> qr(m1); 92 VERIFY(qr.isInjective()); 93 VERIFY(qr.isInvertible()); 94 VERIFY(qr.isSurjective()); 95 96 check_solverbase<MatrixType, MatrixType>(m1, qr, size, size, size); 97 98 // now construct a matrix with prescribed determinant 99 m1.setZero(); 100 for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>(); 101 RealScalar absdet = abs(m1.diagonal().prod()); 102 m3 = qr.matrixQ(); // get a unitary 103 m1 = m3 * m1 * m3; 104 qr.compute(m1); 105 VERIFY_IS_APPROX(absdet, qr.absDeterminant()); 106 VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); 107 } 108 109 template<typename MatrixType> void qr_verify_assert() 110 { 111 MatrixType tmp; 112 113 FullPivHouseholderQR<MatrixType> qr; 114 VERIFY_RAISES_ASSERT(qr.matrixQR()) 115 VERIFY_RAISES_ASSERT(qr.solve(tmp)) 116 VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp)) 117 VERIFY_RAISES_ASSERT(qr.adjoint().solve(tmp)) 118 VERIFY_RAISES_ASSERT(qr.matrixQ()) 119 VERIFY_RAISES_ASSERT(qr.dimensionOfKernel()) 120 VERIFY_RAISES_ASSERT(qr.isInjective()) 121 VERIFY_RAISES_ASSERT(qr.isSurjective()) 122 VERIFY_RAISES_ASSERT(qr.isInvertible()) 123 VERIFY_RAISES_ASSERT(qr.inverse()) 124 VERIFY_RAISES_ASSERT(qr.absDeterminant()) 125 VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) 126 } 127 128 EIGEN_DECLARE_TEST(qr_fullpivoting) 129 { 130 for(int i = 0; i < 1; i++) { 131 CALL_SUBTEST_5( qr<Matrix3f>() ); 132 CALL_SUBTEST_6( qr<Matrix3d>() ); 133 CALL_SUBTEST_8( qr<Matrix2f>() ); 134 CALL_SUBTEST_1( qr<MatrixXf>() ); 135 CALL_SUBTEST_2( qr<MatrixXd>() ); 136 CALL_SUBTEST_3( qr<MatrixXcd>() ); 137 } 138 139 for(int i = 0; i < g_repeat; i++) { 140 CALL_SUBTEST_1( qr_invertible<MatrixXf>() ); 141 CALL_SUBTEST_2( qr_invertible<MatrixXd>() ); 142 CALL_SUBTEST_4( qr_invertible<MatrixXcf>() ); 143 CALL_SUBTEST_3( qr_invertible<MatrixXcd>() ); 144 } 145 146 CALL_SUBTEST_5(qr_verify_assert<Matrix3f>()); 147 CALL_SUBTEST_6(qr_verify_assert<Matrix3d>()); 148 CALL_SUBTEST_1(qr_verify_assert<MatrixXf>()); 149 CALL_SUBTEST_2(qr_verify_assert<MatrixXd>()); 150 CALL_SUBTEST_4(qr_verify_assert<MatrixXcf>()); 151 CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>()); 152 153 // Test problem size constructors 154 CALL_SUBTEST_7(FullPivHouseholderQR<MatrixXf>(10, 20)); 155 CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(10,20))); 156 CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(Matrix<float,10,20>::Random()))); 157 CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(20,10))); 158 CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(Matrix<float,20,10>::Random()))); 159 }