Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame^] | 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 | |
| 14 | template<typename MatrixType> void qr() |
| 15 | { |
| 16 | typedef typename MatrixType::Index Index; |
| 17 | |
| 18 | Index rows = internal::random<Index>(20,200), cols = internal::random<int>(20,200), cols2 = internal::random<int>(20,200); |
| 19 | Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1); |
| 20 | |
| 21 | typedef typename MatrixType::Scalar Scalar; |
| 22 | typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; |
| 23 | typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; |
| 24 | MatrixType m1; |
| 25 | createRandomPIMatrixOfRank(rank,rows,cols,m1); |
| 26 | FullPivHouseholderQR<MatrixType> qr(m1); |
| 27 | VERIFY(rank == qr.rank()); |
| 28 | VERIFY(cols - qr.rank() == qr.dimensionOfKernel()); |
| 29 | VERIFY(!qr.isInjective()); |
| 30 | VERIFY(!qr.isInvertible()); |
| 31 | VERIFY(!qr.isSurjective()); |
| 32 | |
| 33 | MatrixType r = qr.matrixQR(); |
| 34 | |
| 35 | MatrixQType q = qr.matrixQ(); |
| 36 | VERIFY_IS_UNITARY(q); |
| 37 | |
| 38 | // FIXME need better way to construct trapezoid |
| 39 | for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); |
| 40 | |
| 41 | MatrixType c = qr.matrixQ() * r * qr.colsPermutation().inverse(); |
| 42 | |
| 43 | VERIFY_IS_APPROX(m1, c); |
| 44 | |
| 45 | MatrixType m2 = MatrixType::Random(cols,cols2); |
| 46 | MatrixType m3 = m1*m2; |
| 47 | m2 = MatrixType::Random(cols,cols2); |
| 48 | m2 = qr.solve(m3); |
| 49 | VERIFY_IS_APPROX(m3, m1*m2); |
| 50 | } |
| 51 | |
| 52 | template<typename MatrixType> void qr_invertible() |
| 53 | { |
| 54 | typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; |
| 55 | typedef typename MatrixType::Scalar Scalar; |
| 56 | |
| 57 | int size = internal::random<int>(10,50); |
| 58 | |
| 59 | MatrixType m1(size, size), m2(size, size), m3(size, size); |
| 60 | m1 = MatrixType::Random(size,size); |
| 61 | |
| 62 | if (internal::is_same<RealScalar,float>::value) |
| 63 | { |
| 64 | // let's build a matrix more stable to inverse |
| 65 | MatrixType a = MatrixType::Random(size,size*2); |
| 66 | m1 += a * a.adjoint(); |
| 67 | } |
| 68 | |
| 69 | FullPivHouseholderQR<MatrixType> qr(m1); |
| 70 | VERIFY(qr.isInjective()); |
| 71 | VERIFY(qr.isInvertible()); |
| 72 | VERIFY(qr.isSurjective()); |
| 73 | |
| 74 | m3 = MatrixType::Random(size,size); |
| 75 | m2 = qr.solve(m3); |
| 76 | VERIFY_IS_APPROX(m3, m1*m2); |
| 77 | |
| 78 | // now construct a matrix with prescribed determinant |
| 79 | m1.setZero(); |
| 80 | for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>(); |
| 81 | RealScalar absdet = internal::abs(m1.diagonal().prod()); |
| 82 | m3 = qr.matrixQ(); // get a unitary |
| 83 | m1 = m3 * m1 * m3; |
| 84 | qr.compute(m1); |
| 85 | VERIFY_IS_APPROX(absdet, qr.absDeterminant()); |
| 86 | VERIFY_IS_APPROX(internal::log(absdet), qr.logAbsDeterminant()); |
| 87 | } |
| 88 | |
| 89 | template<typename MatrixType> void qr_verify_assert() |
| 90 | { |
| 91 | MatrixType tmp; |
| 92 | |
| 93 | FullPivHouseholderQR<MatrixType> qr; |
| 94 | VERIFY_RAISES_ASSERT(qr.matrixQR()) |
| 95 | VERIFY_RAISES_ASSERT(qr.solve(tmp)) |
| 96 | VERIFY_RAISES_ASSERT(qr.matrixQ()) |
| 97 | VERIFY_RAISES_ASSERT(qr.dimensionOfKernel()) |
| 98 | VERIFY_RAISES_ASSERT(qr.isInjective()) |
| 99 | VERIFY_RAISES_ASSERT(qr.isSurjective()) |
| 100 | VERIFY_RAISES_ASSERT(qr.isInvertible()) |
| 101 | VERIFY_RAISES_ASSERT(qr.inverse()) |
| 102 | VERIFY_RAISES_ASSERT(qr.absDeterminant()) |
| 103 | VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) |
| 104 | } |
| 105 | |
| 106 | void test_qr_fullpivoting() |
| 107 | { |
| 108 | for(int i = 0; i < 1; i++) { |
| 109 | // FIXME : very weird bug here |
| 110 | // CALL_SUBTEST(qr(Matrix2f()) ); |
| 111 | CALL_SUBTEST_1( qr<MatrixXf>() ); |
| 112 | CALL_SUBTEST_2( qr<MatrixXd>() ); |
| 113 | CALL_SUBTEST_3( qr<MatrixXcd>() ); |
| 114 | } |
| 115 | |
| 116 | for(int i = 0; i < g_repeat; i++) { |
| 117 | CALL_SUBTEST_1( qr_invertible<MatrixXf>() ); |
| 118 | CALL_SUBTEST_2( qr_invertible<MatrixXd>() ); |
| 119 | CALL_SUBTEST_4( qr_invertible<MatrixXcf>() ); |
| 120 | CALL_SUBTEST_3( qr_invertible<MatrixXcd>() ); |
| 121 | } |
| 122 | |
| 123 | CALL_SUBTEST_5(qr_verify_assert<Matrix3f>()); |
| 124 | CALL_SUBTEST_6(qr_verify_assert<Matrix3d>()); |
| 125 | CALL_SUBTEST_1(qr_verify_assert<MatrixXf>()); |
| 126 | CALL_SUBTEST_2(qr_verify_assert<MatrixXd>()); |
| 127 | CALL_SUBTEST_4(qr_verify_assert<MatrixXcf>()); |
| 128 | CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>()); |
| 129 | |
| 130 | // Test problem size constructors |
| 131 | CALL_SUBTEST_7(FullPivHouseholderQR<MatrixXf>(10, 20)); |
| 132 | } |