ground/gcs/src/libs/eigen/test/inplace_decomposition.cpp

   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 void 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 }