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 17template<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 82void 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} 111