lu.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-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "main.h"
00026 #include <Eigen/LU>
00027 using namespace std;
00028 
00029 template<typename MatrixType> void lu_non_invertible()
00030 {
00031   typedef typename MatrixType::Index Index;
00032   typedef typename MatrixType::Scalar Scalar;
00033   typedef typename MatrixType::RealScalar RealScalar;
00034   /* this test covers the following files:
00035      LU.h
00036   */
00037   Index rows, cols, cols2;
00038   if(MatrixType::RowsAtCompileTime==Dynamic)
00039   {
00040     rows = internal::random<Index>(2,200);
00041   }
00042   else
00043   {
00044     rows = MatrixType::RowsAtCompileTime;
00045   }
00046   if(MatrixType::ColsAtCompileTime==Dynamic)
00047   {
00048     cols = internal::random<Index>(2,200);
00049     cols2 = internal::random<int>(2,200);
00050   }
00051   else
00052   {
00053     cols2 = cols = MatrixType::ColsAtCompileTime;
00054   }
00055 
00056   enum {
00057     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00058     ColsAtCompileTime = MatrixType::ColsAtCompileTime
00059   };
00060   typedef typename internal::kernel_retval_base<FullPivLU<MatrixType> >::ReturnType KernelMatrixType;
00061   typedef typename internal::image_retval_base<FullPivLU<MatrixType> >::ReturnType ImageMatrixType;
00062   typedef Matrix<typename MatrixType::Scalar, ColsAtCompileTime, ColsAtCompileTime>
00063           CMatrixType;
00064   typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, RowsAtCompileTime>
00065           RMatrixType;
00066 
00067   Index rank = internal::random<Index>(1, std::min(rows, cols)-1);
00068 
00069   // The image of the zero matrix should consist of a single (zero) column vector
00070   VERIFY((MatrixType::Zero(rows,cols).fullPivLu().image(MatrixType::Zero(rows,cols)).cols() == 1));
00071 
00072   MatrixType m1(rows, cols), m3(rows, cols2);
00073   CMatrixType m2(cols, cols2);
00074   createRandomPIMatrixOfRank(rank, rows, cols, m1);
00075 
00076   FullPivLU<MatrixType> lu;
00077 
00078   // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank
00079   // of singular values are either 0 or 1.
00080   // So it's not clear at all that the epsilon should play any role there.
00081   lu.setThreshold(RealScalar(0.01));
00082   lu.compute(m1);
00083 
00084   MatrixType u(rows,cols);
00085   u = lu.matrixLU().template triangularView<Upper>();
00086   RMatrixType l = RMatrixType::Identity(rows,rows);
00087   l.block(0,0,rows,std::min(rows,cols)).template triangularView<StrictlyLower>()
00088     = lu.matrixLU().block(0,0,rows,std::min(rows,cols));
00089 
00090   VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u);
00091 
00092   KernelMatrixType m1kernel = lu.kernel();
00093   ImageMatrixType m1image = lu.image(m1);
00094 
00095   VERIFY_IS_APPROX(m1, lu.reconstructedMatrix());
00096   VERIFY(rank == lu.rank());
00097   VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
00098   VERIFY(!lu.isInjective());
00099   VERIFY(!lu.isInvertible());
00100   VERIFY(!lu.isSurjective());
00101   VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
00102   VERIFY(m1image.fullPivLu().rank() == rank);
00103   VERIFY_IS_APPROX(m1 * m1.adjoint() * m1image, m1image);
00104 
00105   m2 = CMatrixType::Random(cols,cols2);
00106   m3 = m1*m2;
00107   m2 = CMatrixType::Random(cols,cols2);
00108   // test that the code, which does resize(), may be applied to an xpr
00109   m2.block(0,0,m2.rows(),m2.cols()) = lu.solve(m3);
00110   VERIFY_IS_APPROX(m3, m1*m2);
00111 }
00112 
00113 template<typename MatrixType> void lu_invertible()
00114 {
00115   /* this test covers the following files:
00116      LU.h
00117   */
00118   typedef typename MatrixType::Scalar Scalar;
00119   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00120   int size = internal::random<int>(1,200);
00121 
00122   MatrixType m1(size, size), m2(size, size), m3(size, size);
00123   FullPivLU<MatrixType> lu;
00124   lu.setThreshold(RealScalar(0.01));
00125   do {
00126     m1 = MatrixType::Random(size,size);
00127     lu.compute(m1);
00128   } while(!lu.isInvertible());
00129 
00130   VERIFY_IS_APPROX(m1, lu.reconstructedMatrix());
00131   VERIFY(0 == lu.dimensionOfKernel());
00132   VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector
00133   VERIFY(size == lu.rank());
00134   VERIFY(lu.isInjective());
00135   VERIFY(lu.isSurjective());
00136   VERIFY(lu.isInvertible());
00137   VERIFY(lu.image(m1).fullPivLu().isInvertible());
00138   m3 = MatrixType::Random(size,size);
00139   m2 = lu.solve(m3);
00140   VERIFY_IS_APPROX(m3, m1*m2);
00141   VERIFY_IS_APPROX(m2, lu.inverse()*m3);
00142 }
00143 
00144 template<typename MatrixType> void lu_partial_piv()
00145 {
00146   /* this test covers the following files:
00147      PartialPivLU.h
00148   */
00149   typedef typename MatrixType::Index Index;
00150   typedef typename MatrixType::Scalar Scalar;
00151   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00152   Index rows = internal::random<Index>(1,4);
00153   Index cols = rows;
00154 
00155   MatrixType m1(cols, rows);
00156   m1.setRandom();
00157   PartialPivLU<MatrixType> plu(m1);
00158 
00159   VERIFY_IS_APPROX(m1, plu.reconstructedMatrix());
00160 }
00161 
00162 template<typename MatrixType> void lu_verify_assert()
00163 {
00164   MatrixType tmp;
00165 
00166   FullPivLU<MatrixType> lu;
00167   VERIFY_RAISES_ASSERT(lu.matrixLU())
00168   VERIFY_RAISES_ASSERT(lu.permutationP())
00169   VERIFY_RAISES_ASSERT(lu.permutationQ())
00170   VERIFY_RAISES_ASSERT(lu.kernel())
00171   VERIFY_RAISES_ASSERT(lu.image(tmp))
00172   VERIFY_RAISES_ASSERT(lu.solve(tmp))
00173   VERIFY_RAISES_ASSERT(lu.determinant())
00174   VERIFY_RAISES_ASSERT(lu.rank())
00175   VERIFY_RAISES_ASSERT(lu.dimensionOfKernel())
00176   VERIFY_RAISES_ASSERT(lu.isInjective())
00177   VERIFY_RAISES_ASSERT(lu.isSurjective())
00178   VERIFY_RAISES_ASSERT(lu.isInvertible())
00179   VERIFY_RAISES_ASSERT(lu.inverse())
00180 
00181   PartialPivLU<MatrixType> plu;
00182   VERIFY_RAISES_ASSERT(plu.matrixLU())
00183   VERIFY_RAISES_ASSERT(plu.permutationP())
00184   VERIFY_RAISES_ASSERT(plu.solve(tmp))
00185   VERIFY_RAISES_ASSERT(plu.determinant())
00186   VERIFY_RAISES_ASSERT(plu.inverse())
00187 }
00188 
00189 void test_lu()
00190 {
00191   for(int i = 0; i < g_repeat; i++) {
00192     CALL_SUBTEST_1( lu_non_invertible<Matrix3f>() );
00193     CALL_SUBTEST_1( lu_verify_assert<Matrix3f>() );
00194 
00195     CALL_SUBTEST_2( (lu_non_invertible<Matrix<double, 4, 6> >()) );
00196     CALL_SUBTEST_2( (lu_verify_assert<Matrix<double, 4, 6> >()) );
00197 
00198     CALL_SUBTEST_3( lu_non_invertible<MatrixXf>() );
00199     CALL_SUBTEST_3( lu_invertible<MatrixXf>() );
00200     CALL_SUBTEST_3( lu_verify_assert<MatrixXf>() );
00201 
00202     CALL_SUBTEST_4( lu_non_invertible<MatrixXd>() );
00203     CALL_SUBTEST_4( lu_invertible<MatrixXd>() );
00204     CALL_SUBTEST_4( lu_partial_piv<MatrixXd>() );
00205     CALL_SUBTEST_4( lu_verify_assert<MatrixXd>() );
00206 
00207     CALL_SUBTEST_5( lu_non_invertible<MatrixXcf>() );
00208     CALL_SUBTEST_5( lu_invertible<MatrixXcf>() );
00209     CALL_SUBTEST_5( lu_verify_assert<MatrixXcf>() );
00210 
00211     CALL_SUBTEST_6( lu_non_invertible<MatrixXcd>() );
00212     CALL_SUBTEST_6( lu_invertible<MatrixXcd>() );
00213     CALL_SUBTEST_6( lu_partial_piv<MatrixXcd>() );
00214     CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() );
00215 
00216     CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() ));
00217 
00218     // Test problem size constructors
00219     CALL_SUBTEST_9( PartialPivLU<MatrixXf>(10) );
00220     CALL_SUBTEST_9( FullPivLU<MatrixXf>(10, 20); );
00221   }
00222 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:58