blob: 92d0d91b6b42e0d27031f968af38e29233df8e23 [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2016 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #include "main.h" #include #include #include // This file test inplace decomposition through Ref<>, as supported by Cholesky, LU, and QR decompositions. template void inplace(bool square = false, bool SPD = false) { typedef typename MatrixType::Scalar Scalar; typedef Matrix RhsType; typedef Matrix ResType; Index rows = MatrixType::RowsAtCompileTime==Dynamic ? internal::random(2,EIGEN_TEST_MAX_SIZE/2) : Index(MatrixType::RowsAtCompileTime); Index cols = MatrixType::ColsAtCompileTime==Dynamic ? (square?rows:internal::random(2,rows)) : Index(MatrixType::ColsAtCompileTime); MatrixType A = MatrixType::Random(rows,cols); RhsType b = RhsType::Random(rows); ResType x(cols); if(SPD) { assert(square); A.topRows(cols) = A.topRows(cols).adjoint() * A.topRows(cols); A.diagonal().array() += 1e-3; } MatrixType A0 = A; MatrixType A1 = A; DecType dec(A); // Check that the content of A has been modified VERIFY_IS_NOT_APPROX( A, A0 ); // Check that the decomposition is correct: if(rows==cols) { VERIFY_IS_APPROX( A0 * (x = dec.solve(b)), b ); } else { VERIFY_IS_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); } // Check that modifying A breaks the current dec: A.setRandom(); if(rows==cols) { VERIFY_IS_NOT_APPROX( A0 * (x = dec.solve(b)), b ); } else { VERIFY_IS_NOT_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); } // Check that calling compute(A1) does not modify A1: A = A0; dec.compute(A1); VERIFY_IS_EQUAL(A0,A1); VERIFY_IS_NOT_APPROX( A, A0 ); if(rows==cols) { VERIFY_IS_APPROX( A0 * (x = dec.solve(b)), b ); } else { VERIFY_IS_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); } } void test_inplace_decomposition() { EIGEN_UNUSED typedef Matrix Matrix43d; for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(( inplace >, MatrixXd>(true,true) )); CALL_SUBTEST_1(( inplace >, Matrix4d>(true,true) )); CALL_SUBTEST_2(( inplace >, MatrixXd>(true,true) )); CALL_SUBTEST_2(( inplace >, Matrix4d>(true,true) )); CALL_SUBTEST_3(( inplace >, MatrixXd>(true,false) )); CALL_SUBTEST_3(( inplace >, Matrix4d>(true,false) )); CALL_SUBTEST_4(( inplace >, MatrixXd>(true,false) )); CALL_SUBTEST_4(( inplace >, Matrix4d>(true,false) )); CALL_SUBTEST_5(( inplace >, MatrixXd>(false,false) )); CALL_SUBTEST_5(( inplace >, Matrix43d>(false,false) )); CALL_SUBTEST_6(( inplace >, MatrixXd>(false,false) )); CALL_SUBTEST_6(( inplace >, Matrix43d>(false,false) )); CALL_SUBTEST_7(( inplace >, MatrixXd>(false,false) )); CALL_SUBTEST_7(( inplace >, Matrix43d>(false,false) )); CALL_SUBTEST_8(( inplace >, MatrixXd>(false,false) )); CALL_SUBTEST_8(( inplace >, Matrix43d>(false,false) )); } }