00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "main.h"
00027 #include <Eigen/QR>
00028
00029 template<typename MatrixType> void qr()
00030 {
00031 typedef typename MatrixType::Index Index;
00032
00033 Index rows = internal::random<Index>(2,200), cols = internal::random<Index>(2,200), cols2 = internal::random<Index>(2,200);
00034 Index rank = internal::random<Index>(1, std::min(rows, cols)-1);
00035
00036 typedef typename MatrixType::Scalar Scalar;
00037 typedef typename MatrixType::RealScalar RealScalar;
00038 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
00039 typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
00040 MatrixType m1;
00041 createRandomPIMatrixOfRank(rank,rows,cols,m1);
00042 ColPivHouseholderQR<MatrixType> qr(m1);
00043 VERIFY(rank == qr.rank());
00044 VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
00045 VERIFY(!qr.isInjective());
00046 VERIFY(!qr.isInvertible());
00047 VERIFY(!qr.isSurjective());
00048
00049 MatrixQType q = qr.householderQ();
00050 VERIFY_IS_UNITARY(q);
00051
00052 MatrixType r = qr.matrixQR().template triangularView<Upper>();
00053 MatrixType c = q * r * qr.colsPermutation().inverse();
00054 VERIFY_IS_APPROX(m1, c);
00055
00056 MatrixType m2 = MatrixType::Random(cols,cols2);
00057 MatrixType m3 = m1*m2;
00058 m2 = MatrixType::Random(cols,cols2);
00059 m2 = qr.solve(m3);
00060 VERIFY_IS_APPROX(m3, m1*m2);
00061 }
00062
00063 template<typename MatrixType, int Cols2> void qr_fixedsize()
00064 {
00065 enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
00066 typedef typename MatrixType::Scalar Scalar;
00067 int rank = internal::random<int>(1, std::min(int(Rows), int(Cols))-1);
00068 Matrix<Scalar,Rows,Cols> m1;
00069 createRandomPIMatrixOfRank(rank,Rows,Cols,m1);
00070 ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
00071 VERIFY(rank == qr.rank());
00072 VERIFY(Cols - qr.rank() == qr.dimensionOfKernel());
00073 VERIFY(qr.isInjective() == (rank == Rows));
00074 VERIFY(qr.isSurjective() == (rank == Cols));
00075 VERIFY(qr.isInvertible() == (qr.isInjective() && qr.isSurjective()));
00076
00077 Matrix<Scalar,Rows,Cols> r = qr.matrixQR().template triangularView<Upper>();
00078 Matrix<Scalar,Rows,Cols> c = qr.householderQ() * r * qr.colsPermutation().inverse();
00079 VERIFY_IS_APPROX(m1, c);
00080
00081 Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
00082 Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
00083 m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
00084 m2 = qr.solve(m3);
00085 VERIFY_IS_APPROX(m3, m1*m2);
00086 }
00087
00088 template<typename MatrixType> void qr_invertible()
00089 {
00090 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00091 typedef typename MatrixType::Scalar Scalar;
00092
00093 int size = internal::random<int>(10,50);
00094
00095 MatrixType m1(size, size), m2(size, size), m3(size, size);
00096 m1 = MatrixType::Random(size,size);
00097
00098 if (internal::is_same<RealScalar,float>::value)
00099 {
00100
00101 MatrixType a = MatrixType::Random(size,size*2);
00102 m1 += a * a.adjoint();
00103 }
00104
00105 ColPivHouseholderQR<MatrixType> qr(m1);
00106 m3 = MatrixType::Random(size,size);
00107 m2 = qr.solve(m3);
00108
00109
00110
00111 m1.setZero();
00112 for(int i = 0; i < size; i++) m1(i,i) = internal::random<Scalar>();
00113 RealScalar absdet = internal::abs(m1.diagonal().prod());
00114 m3 = qr.householderQ();
00115 m1 = m3 * m1 * m3;
00116 qr.compute(m1);
00117 VERIFY_IS_APPROX(absdet, qr.absDeterminant());
00118 VERIFY_IS_APPROX(internal::log(absdet), qr.logAbsDeterminant());
00119 }
00120
00121 template<typename MatrixType> void qr_verify_assert()
00122 {
00123 MatrixType tmp;
00124
00125 ColPivHouseholderQR<MatrixType> qr;
00126 VERIFY_RAISES_ASSERT(qr.matrixQR())
00127 VERIFY_RAISES_ASSERT(qr.solve(tmp))
00128 VERIFY_RAISES_ASSERT(qr.householderQ())
00129 VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
00130 VERIFY_RAISES_ASSERT(qr.isInjective())
00131 VERIFY_RAISES_ASSERT(qr.isSurjective())
00132 VERIFY_RAISES_ASSERT(qr.isInvertible())
00133 VERIFY_RAISES_ASSERT(qr.inverse())
00134 VERIFY_RAISES_ASSERT(qr.absDeterminant())
00135 VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
00136 }
00137
00138 void test_qr_colpivoting()
00139 {
00140 for(int i = 0; i < g_repeat; i++) {
00141 CALL_SUBTEST_1( qr<MatrixXf>() );
00142 CALL_SUBTEST_2( qr<MatrixXd>() );
00143 CALL_SUBTEST_3( qr<MatrixXcd>() );
00144 CALL_SUBTEST_4(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
00145 CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
00146 CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,1,1>, 1 >() ));
00147 }
00148
00149 for(int i = 0; i < g_repeat; i++) {
00150 CALL_SUBTEST_1( qr_invertible<MatrixXf>() );
00151 CALL_SUBTEST_2( qr_invertible<MatrixXd>() );
00152 CALL_SUBTEST_6( qr_invertible<MatrixXcf>() );
00153 CALL_SUBTEST_3( qr_invertible<MatrixXcd>() );
00154 }
00155
00156 CALL_SUBTEST_7(qr_verify_assert<Matrix3f>());
00157 CALL_SUBTEST_8(qr_verify_assert<Matrix3d>());
00158 CALL_SUBTEST_1(qr_verify_assert<MatrixXf>());
00159 CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
00160 CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>());
00161 CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
00162
00163
00164 CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20));
00165 }