cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

inplace_decomposition.cpp (3880B)

      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #include "main.h"
     11 #include <Eigen/LU>
     12 #include <Eigen/Cholesky>
     13 #include <Eigen/QR>
     14 
     15 // This file test inplace decomposition through Ref<>, as supported by Cholesky, LU, and QR decompositions.
     16 
     17 template<typename DecType,typename MatrixType> void inplace(bool square = false, bool SPD = false)
     18 {
     19   typedef typename MatrixType::Scalar Scalar;
     20   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> RhsType;
     21   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> ResType;
     22 
     23   Index rows = MatrixType::RowsAtCompileTime==Dynamic ? internal::random<Index>(2,EIGEN_TEST_MAX_SIZE/2) : Index(MatrixType::RowsAtCompileTime);
     24   Index cols = MatrixType::ColsAtCompileTime==Dynamic ? (square?rows:internal::random<Index>(2,rows))    : Index(MatrixType::ColsAtCompileTime);
     25 
     26   MatrixType A = MatrixType::Random(rows,cols);
     27   RhsType b = RhsType::Random(rows);
     28   ResType x(cols);
     29 
     30   if(SPD)
     31   {
     32     assert(square);
     33     A.topRows(cols) = A.topRows(cols).adjoint() * A.topRows(cols);
     34     A.diagonal().array() += 1e-3;
     35   }
     36 
     37   MatrixType A0 = A;
     38   MatrixType A1 = A;
     39 
     40   DecType dec(A);
     41 
     42   // Check that the content of A has been modified
     43   VERIFY_IS_NOT_APPROX( A, A0 );
     44 
     45   // Check that the decomposition is correct:
     46   if(rows==cols)
     47   {
     48     VERIFY_IS_APPROX( A0 * (x = dec.solve(b)), b );
     49   }
     50   else
     51   {
     52     VERIFY_IS_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b );
     53   }
     54 
     55   // Check that modifying A breaks the current dec:
     56   A.setRandom();
     57   if(rows==cols)
     58   {
     59     VERIFY_IS_NOT_APPROX( A0 * (x = dec.solve(b)), b );
     60   }
     61   else
     62   {
     63     VERIFY_IS_NOT_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b );
     64   }
     65 
     66   // Check that calling compute(A1) does not modify A1:
     67   A = A0;
     68   dec.compute(A1);
     69   VERIFY_IS_EQUAL(A0,A1);
     70   VERIFY_IS_NOT_APPROX( A, A0 );
     71   if(rows==cols)
     72   {
     73     VERIFY_IS_APPROX( A0 * (x = dec.solve(b)), b );
     74   }
     75   else
     76   {
     77     VERIFY_IS_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b );
     78   }
     79 }
     80 
     81 
     82 EIGEN_DECLARE_TEST(inplace_decomposition)
     83 {
     84   EIGEN_UNUSED typedef Matrix<double,4,3> Matrix43d;
     85   for(int i = 0; i < g_repeat; i++) {
     86     CALL_SUBTEST_1(( inplace<LLT<Ref<MatrixXd> >, MatrixXd>(true,true) ));
     87     CALL_SUBTEST_1(( inplace<LLT<Ref<Matrix4d> >, Matrix4d>(true,true) ));
     88 
     89     CALL_SUBTEST_2(( inplace<LDLT<Ref<MatrixXd> >, MatrixXd>(true,true) ));
     90     CALL_SUBTEST_2(( inplace<LDLT<Ref<Matrix4d> >, Matrix4d>(true,true) ));
     91 
     92     CALL_SUBTEST_3(( inplace<PartialPivLU<Ref<MatrixXd> >, MatrixXd>(true,false) ));
     93     CALL_SUBTEST_3(( inplace<PartialPivLU<Ref<Matrix4d> >, Matrix4d>(true,false) ));
     94 
     95     CALL_SUBTEST_4(( inplace<FullPivLU<Ref<MatrixXd> >, MatrixXd>(true,false) ));
     96     CALL_SUBTEST_4(( inplace<FullPivLU<Ref<Matrix4d> >, Matrix4d>(true,false) ));
     97 
     98     CALL_SUBTEST_5(( inplace<HouseholderQR<Ref<MatrixXd> >, MatrixXd>(false,false) ));
     99     CALL_SUBTEST_5(( inplace<HouseholderQR<Ref<Matrix43d> >, Matrix43d>(false,false) ));
    100 
    101     CALL_SUBTEST_6(( inplace<ColPivHouseholderQR<Ref<MatrixXd> >, MatrixXd>(false,false) ));
    102     CALL_SUBTEST_6(( inplace<ColPivHouseholderQR<Ref<Matrix43d> >, Matrix43d>(false,false) ));
    103 
    104     CALL_SUBTEST_7(( inplace<FullPivHouseholderQR<Ref<MatrixXd> >, MatrixXd>(false,false) ));
    105     CALL_SUBTEST_7(( inplace<FullPivHouseholderQR<Ref<Matrix43d> >, Matrix43d>(false,false) ));
    106 
    107     CALL_SUBTEST_8(( inplace<CompleteOrthogonalDecomposition<Ref<MatrixXd> >, MatrixXd>(false,false) ));
    108     CALL_SUBTEST_8(( inplace<CompleteOrthogonalDecomposition<Ref<Matrix43d> >, Matrix43d>(false,false) ));
    109   }
    110 }