qr_colpivoting.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
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     // let's build a matrix more stable to inverse
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   //VERIFY_IS_APPROX(m3, m1*m2);
00109 
00110   // now construct a matrix with prescribed determinant
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(); // get a unitary
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   // Test problem size constructors
00164   CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20));
00165 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:15