inplace_decomposition.cpp
Go to the documentation of this file.
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;
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 
109  }
110 }
static Matrix A1
Robust Cholesky decomposition of a matrix with pivoting.
Definition: LDLT.h:59
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define CALL_SUBTEST_6(FUNC)
#define CALL_SUBTEST_4(FUNC)
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Scalar * b
Definition: benchVecAdd.cpp:17
#define CALL_SUBTEST_3(FUNC)
#define CALL_SUBTEST_7(FUNC)
LU decomposition of a matrix with partial pivoting, and related features.
MatrixXf MatrixType
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
EIGEN_DECLARE_TEST(inplace_decomposition)
Complete orthogonal decomposition (COD) of a matrix.
#define EIGEN_UNUSED
Definition: Macros.h:1067
#define VERIFY_IS_APPROX(a, b)
EIGEN_DEVICE_FUNC const SquareReturnType square() const
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
#define CALL_SUBTEST_1(FUNC)
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
Standard Cholesky decomposition (LL^T) of a matrix and associated features.
Definition: LLT.h:66
static int g_repeat
Definition: main.h:169
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
#define CALL_SUBTEST_8(FUNC)
Array< double, 1, 3 > e(1./3., 0.5, 2.)
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
#define CALL_SUBTEST_5(FUNC)
#define EIGEN_TEST_MAX_SIZE
LU decomposition of a matrix with complete pivoting, and related features.
Householder QR decomposition of a matrix.
#define VERIFY_IS_NOT_APPROX(a, b)
#define CALL_SUBTEST_2(FUNC)
const int Dynamic
Definition: Constants.h:22
void inplace(bool square=false, bool SPD=false)
The matrix class, also used for vectors and row-vectors.
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:22