cholesky.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 //
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 #ifndef EIGEN_NO_ASSERTION_CHECKING
00026 #define EIGEN_NO_ASSERTION_CHECKING
00027 #endif
00028 
00029 static int nb_temporaries;
00030 
00031 #define EIGEN_DENSE_STORAGE_CTOR_PLUGIN { if(size!=0) nb_temporaries++; }
00032 
00033 #include "main.h"
00034 #include <Eigen/Cholesky>
00035 #include <Eigen/QR>
00036 
00037 #define VERIFY_EVALUATION_COUNT(XPR,N) {\
00038     nb_temporaries = 0; \
00039     XPR; \
00040     if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \
00041     VERIFY( (#XPR) && nb_temporaries==N ); \
00042   }
00043 
00044 #ifdef HAS_GSL
00045 #include "gsl_helper.h"
00046 #endif
00047 
00048 template<typename MatrixType> void cholesky(const MatrixType& m)
00049 {
00050   typedef typename MatrixType::Index Index;
00051   /* this test covers the following files:
00052      LLT.h LDLT.h
00053   */
00054   Index rows = m.rows();
00055   Index cols = m.cols();
00056 
00057   typedef typename MatrixType::Scalar Scalar;
00058   typedef typename NumTraits<Scalar>::Real RealScalar;
00059   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
00060   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00061 
00062   MatrixType a0 = MatrixType::Random(rows,cols);
00063   VectorType vecB = VectorType::Random(rows), vecX(rows);
00064   MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
00065   SquareMatrixType symm =  a0 * a0.adjoint();
00066   // let's make sure the matrix is not singular or near singular
00067   for (int k=0; k<3; ++k)
00068   {
00069     MatrixType a1 = MatrixType::Random(rows,cols);
00070     symm += a1 * a1.adjoint();
00071   }
00072 
00073   SquareMatrixType symmUp = symm.template triangularView<Upper>();
00074   SquareMatrixType symmLo = symm.template triangularView<Lower>();
00075 
00076   // to test if really Cholesky only uses the upper triangular part, uncomment the following
00077   // FIXME: currently that fails !!
00078   //symm.template part<StrictlyLower>().setZero();
00079 
00080   #ifdef HAS_GSL
00081 //   if (internal::is_same<RealScalar,double>::value)
00082 //   {
00083 //     typedef GslTraits<Scalar> Gsl;
00084 //     typename Gsl::Matrix gMatA=0, gSymm=0;
00085 //     typename Gsl::Vector gVecB=0, gVecX=0;
00086 //     convert<MatrixType>(symm, gSymm);
00087 //     convert<MatrixType>(symm, gMatA);
00088 //     convert<VectorType>(vecB, gVecB);
00089 //     convert<VectorType>(vecB, gVecX);
00090 //     Gsl::cholesky(gMatA);
00091 //     Gsl::cholesky_solve(gMatA, gVecB, gVecX);
00092 //     VectorType vecX(rows), _vecX, _vecB;
00093 //     convert(gVecX, _vecX);
00094 //     symm.llt().solve(vecB, &vecX);
00095 //     Gsl::prod(gSymm, gVecX, gVecB);
00096 //     convert(gVecB, _vecB);
00097 //     // test gsl itself !
00098 //     VERIFY_IS_APPROX(vecB, _vecB);
00099 //     VERIFY_IS_APPROX(vecX, _vecX);
00100 //
00101 //     Gsl::free(gMatA);
00102 //     Gsl::free(gSymm);
00103 //     Gsl::free(gVecB);
00104 //     Gsl::free(gVecX);
00105 //   }
00106   #endif
00107 
00108   {
00109     LLT<SquareMatrixType,Lower> chollo(symmLo);
00110     VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
00111     vecX = chollo.solve(vecB);
00112     VERIFY_IS_APPROX(symm * vecX, vecB);
00113     matX = chollo.solve(matB);
00114     VERIFY_IS_APPROX(symm * matX, matB);
00115 
00116     // test the upper mode
00117     LLT<SquareMatrixType,Upper> cholup(symmUp);
00118     VERIFY_IS_APPROX(symm, cholup.reconstructedMatrix());
00119     vecX = cholup.solve(vecB);
00120     VERIFY_IS_APPROX(symm * vecX, vecB);
00121     matX = cholup.solve(matB);
00122     VERIFY_IS_APPROX(symm * matX, matB);
00123 
00124     MatrixType neg = -symmLo;
00125     chollo.compute(neg);
00126     VERIFY(chollo.info()==NumericalIssue);
00127   }
00128 
00129   // LDLT
00130   {
00131     int sign = internal::random<int>()%2 ? 1 : -1;
00132 
00133     if(sign == -1)
00134     {
00135       symm = -symm; // test a negative matrix
00136     }
00137 
00138     SquareMatrixType symmUp = symm.template triangularView<Upper>();
00139     SquareMatrixType symmLo = symm.template triangularView<Lower>();
00140 
00141     LDLT<SquareMatrixType,Lower> ldltlo(symmLo);
00142     VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
00143     vecX = ldltlo.solve(vecB);
00144     VERIFY_IS_APPROX(symm * vecX, vecB);
00145     matX = ldltlo.solve(matB);
00146     VERIFY_IS_APPROX(symm * matX, matB);
00147 
00148     LDLT<SquareMatrixType,Upper> ldltup(symmUp);
00149     VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix());
00150     vecX = ldltup.solve(vecB);
00151     VERIFY_IS_APPROX(symm * vecX, vecB);
00152     matX = ldltup.solve(matB);
00153     VERIFY_IS_APPROX(symm * matX, matB);
00154 
00155     if(MatrixType::RowsAtCompileTime==Dynamic)
00156     {
00157       // note : each inplace permutation requires a small temporary vector (mask)
00158 
00159       // check inplace solve
00160       matX = matB;
00161       VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0);
00162       VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval());
00163 
00164 
00165       matX = matB;
00166       VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
00167       VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
00168     }
00169   }
00170 
00171   // test some special use cases of SelfCwiseBinaryOp:
00172   MatrixType m1 = MatrixType::Random(rows,cols), m2(rows,cols);
00173   m2 = m1;
00174   m2 += symmLo.template selfadjointView<Lower>().llt().solve(matB);
00175   VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
00176   m2 = m1;
00177   m2 -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
00178   VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
00179   m2 = m1;
00180   m2.noalias() += symmLo.template selfadjointView<Lower>().llt().solve(matB);
00181   VERIFY_IS_APPROX(m2, m1 + symmLo.template selfadjointView<Lower>().llt().solve(matB));
00182   m2 = m1;
00183   m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
00184   VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
00185   
00186 }
00187 
00188 template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
00189 {
00190   // classic test
00191   cholesky(m);
00192 
00193   // test mixing real/scalar types
00194 
00195   typedef typename MatrixType::Index Index;
00196 
00197   Index rows = m.rows();
00198   Index cols = m.cols();
00199 
00200   typedef typename MatrixType::Scalar Scalar;
00201   typedef typename NumTraits<Scalar>::Real RealScalar;
00202   typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RealMatrixType;
00203   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00204 
00205   RealMatrixType a0 = RealMatrixType::Random(rows,cols);
00206   VectorType vecB = VectorType::Random(rows), vecX(rows);
00207   MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
00208   RealMatrixType symm =  a0 * a0.adjoint();
00209   // let's make sure the matrix is not singular or near singular
00210   for (int k=0; k<3; ++k)
00211   {
00212     RealMatrixType a1 = RealMatrixType::Random(rows,cols);
00213     symm += a1 * a1.adjoint();
00214   }
00215 
00216   {
00217     RealMatrixType symmLo = symm.template triangularView<Lower>();
00218 
00219     LLT<RealMatrixType,Lower> chollo(symmLo);
00220     VERIFY_IS_APPROX(symm, chollo.reconstructedMatrix());
00221     vecX = chollo.solve(vecB);
00222     VERIFY_IS_APPROX(symm * vecX, vecB);
00223 //     matX = chollo.solve(matB);
00224 //     VERIFY_IS_APPROX(symm * matX, matB);
00225   }
00226 
00227   // LDLT
00228   {
00229     int sign = internal::random<int>()%2 ? 1 : -1;
00230 
00231     if(sign == -1)
00232     {
00233       symm = -symm; // test a negative matrix
00234     }
00235 
00236     RealMatrixType symmLo = symm.template triangularView<Lower>();
00237 
00238     LDLT<RealMatrixType,Lower> ldltlo(symmLo);
00239     VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
00240     vecX = ldltlo.solve(vecB);
00241     VERIFY_IS_APPROX(symm * vecX, vecB);
00242 //     matX = ldltlo.solve(matB);
00243 //     VERIFY_IS_APPROX(symm * matX, matB);
00244   }
00245 
00246 }
00247 
00248 template<typename MatrixType> void cholesky_verify_assert()
00249 {
00250   MatrixType tmp;
00251 
00252   LLT<MatrixType> llt;
00253   VERIFY_RAISES_ASSERT(llt.matrixL())
00254   VERIFY_RAISES_ASSERT(llt.matrixU())
00255   VERIFY_RAISES_ASSERT(llt.solve(tmp))
00256   VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp))
00257 
00258   LDLT<MatrixType> ldlt;
00259   VERIFY_RAISES_ASSERT(ldlt.matrixL())
00260   VERIFY_RAISES_ASSERT(ldlt.permutationP())
00261   VERIFY_RAISES_ASSERT(ldlt.vectorD())
00262   VERIFY_RAISES_ASSERT(ldlt.isPositive())
00263   VERIFY_RAISES_ASSERT(ldlt.isNegative())
00264   VERIFY_RAISES_ASSERT(ldlt.solve(tmp))
00265   VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp))
00266 }
00267 
00268 void test_cholesky()
00269 {
00270   int s;
00271   for(int i = 0; i < g_repeat; i++) {
00272     CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) );
00273     CALL_SUBTEST_3( cholesky(Matrix2d()) );
00274     CALL_SUBTEST_4( cholesky(Matrix3f()) );
00275     CALL_SUBTEST_5( cholesky(Matrix4d()) );
00276     s = internal::random<int>(1,200);
00277     CALL_SUBTEST_2( cholesky(MatrixXd(s,s)) );
00278     s = internal::random<int>(1,100);
00279     CALL_SUBTEST_6( cholesky_cplx(MatrixXcd(s,s)) );
00280   }
00281 
00282   CALL_SUBTEST_4( cholesky_verify_assert<Matrix3f>() );
00283   CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );
00284   CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() );
00285   CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() );
00286 
00287   // Test problem size constructors
00288   CALL_SUBTEST_9( LLT<MatrixXf>(10) );
00289   CALL_SUBTEST_9( LDLT<MatrixXf>(10) );
00290 }


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