From cd3c2451b6915dbd9e27930b1dcb4550a73d8e2e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 11 Oct 2011 13:45:27 +0200 Subject: add a unit test for permutation applied to sparse objects --- test/CMakeLists.txt | 1 + test/main.h | 17 ++++++ test/permutationmatrices.cpp | 17 ------ test/sparse_permutations.cpp | 137 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 155 insertions(+), 17 deletions(-) create mode 100644 test/sparse_permutations.cpp (limited to 'test') diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9fffd02a7..6b0fa1a8e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -129,6 +129,7 @@ ei_add_test(householder) ei_add_test(swap) ei_add_test(conservative_resize) ei_add_test(permutationmatrices) +ei_add_test(sparse_permutations) ei_add_test(eigen2support) ei_add_test(nullary) ei_add_test(nesting_ops "${CMAKE_CXX_FLAGS_DEBUG}") diff --git a/test/main.h b/test/main.h index 2068713b3..d9cc95614 100644 --- a/test/main.h +++ b/test/main.h @@ -355,6 +355,23 @@ void createRandomPIMatrixOfRank(typename MatrixType::Index desired_rank, typenam m = qra.householderQ() * d * qrb.householderQ(); } +template +void randomPermutationVector(PermutationVectorType& v, typename PermutationVectorType::Index size) +{ + typedef typename PermutationVectorType::Index Index; + typedef typename PermutationVectorType::Scalar Scalar; + v.resize(size); + for(Index i = 0; i < size; ++i) v(i) = Scalar(i); + if(size == 1) return; + for(Index n = 0; n < 3 * size; ++n) + { + Index i = internal::random(0, size-1); + Index j; + do j = internal::random(0, size-1); while(j==i); + std::swap(v(i), v(j)); + } +} + } // end namespace Eigen template struct GetDifferentType; diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp index d0fa01310..308838c70 100644 --- a/test/permutationmatrices.cpp +++ b/test/permutationmatrices.cpp @@ -24,23 +24,6 @@ #include "main.h" -template -void randomPermutationVector(PermutationVectorType& v, typename PermutationVectorType::Index size) -{ - typedef typename PermutationVectorType::Index Index; - typedef typename PermutationVectorType::Scalar Scalar; - v.resize(size); - for(Index i = 0; i < size; ++i) v(i) = Scalar(i); - if(size == 1) return; - for(Index n = 0; n < 3 * size; ++n) - { - Index i = internal::random(0, size-1); - Index j; - do j = internal::random(0, size-1); while(j==i); - std::swap(v(i), v(j)); - } -} - using namespace std; template void permutationmatrices(const MatrixType& m) { diff --git a/test/sparse_permutations.cpp b/test/sparse_permutations.cpp new file mode 100644 index 000000000..7f47611a5 --- /dev/null +++ b/test/sparse_permutations.cpp @@ -0,0 +1,137 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#include "sparse.h" + +template void sparse_permutations(const SparseMatrixType& ref) +{ + typedef typename SparseMatrixType::Index Index; + + const Index rows = ref.rows(); + const Index cols = ref.cols(); + typedef typename SparseMatrixType::Scalar Scalar; + typedef Matrix DenseMatrix; + typedef Matrix VectorI; + + double density = (std::max)(8./(rows*cols), 0.01); + + SparseMatrixType mat(rows, cols), up(rows,cols), lo(rows,cols), res; + DenseMatrix mat_d = DenseMatrix::Zero(rows, cols), up_sym_d, lo_sym_d, res_d; + + initSparse(density, mat_d, mat, 0); + + up = mat.template triangularView(); + lo = mat.template triangularView(); + + up_sym_d = mat_d.template selfadjointView(); + lo_sym_d = mat_d.template selfadjointView(); + + VERIFY_IS_APPROX(mat, mat_d); + VERIFY_IS_APPROX(up, DenseMatrix(mat_d.template triangularView())); + VERIFY_IS_APPROX(lo, DenseMatrix(mat_d.template triangularView())); + + PermutationMatrix p, p_null; + VectorI pi; + randomPermutationVector(pi, cols); + p.indices() = pi; + + + res = mat.template selfadjointView().twistedBy(p_null); + res_d = up_sym_d; + VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full"); + + res = mat.template selfadjointView().twistedBy(p_null); + res_d = lo_sym_d; + VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full"); + + + res = up.template selfadjointView().twistedBy(p_null); + res_d = up_sym_d; + VERIFY(res.isApprox(res_d) && "upper selfadjoint to full"); + + res = lo.template selfadjointView().twistedBy(p_null); + res_d = lo_sym_d; + VERIFY(res.isApprox(res_d) && "lower selfadjoint full"); + + + res.template selfadjointView() = mat.template selfadjointView().twistedBy(p); + res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView(); + VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to upper"); + + res.template selfadjointView() = mat.template selfadjointView().twistedBy(p); + res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView(); + VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to upper"); + + res.template selfadjointView() = mat.template selfadjointView().twistedBy(p); + res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView(); + VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to lower"); + + res.template selfadjointView() = mat.template selfadjointView().twistedBy(p); + res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView(); + VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to lower"); + + + res.template selfadjointView() = up.template selfadjointView().twistedBy(p); + res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView(); + VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to upper"); + + res.template selfadjointView() = lo.template selfadjointView().twistedBy(p); + res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView(); + VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to upper"); + + res.template selfadjointView() = lo.template selfadjointView().twistedBy(p); + res_d = ((p * lo_sym_d) * p.inverse()).eval().template triangularView(); + VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to lower"); + + res.template selfadjointView() = up.template selfadjointView().twistedBy(p); + res_d = ((p * up_sym_d) * p.inverse()).eval().template triangularView(); + VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to lower"); + + + res = mat.template selfadjointView().twistedBy(p); + res_d = (p * up_sym_d) * p.inverse(); + VERIFY(res.isApprox(res_d) && "full selfadjoint upper twisted to full"); + + res = mat.template selfadjointView().twistedBy(p); + res_d = (p * lo_sym_d) * p.inverse(); + VERIFY(res.isApprox(res_d) && "full selfadjoint lower twisted to full"); + + res = up.template selfadjointView().twistedBy(p); + res_d = (p * up_sym_d) * p.inverse(); + VERIFY(res.isApprox(res_d) && "upper selfadjoint twisted to full"); + + res = lo.template selfadjointView().twistedBy(p); + res_d = (p * lo_sym_d) * p.inverse(); + VERIFY(res.isApprox(res_d) && "lower selfadjoint twisted to full"); +} + +void test_sparse_permutations() +{ + for(int i = 0; i < g_repeat; i++) { + int s = Eigen::internal::random(1,50); + CALL_SUBTEST_1(( sparse_permutations(SparseMatrix(8, 8)) )); + CALL_SUBTEST_2(( sparse_permutations(SparseMatrix >(s, s)) )); + CALL_SUBTEST_1(( sparse_permutations(SparseMatrix(s, s)) )); + } +} -- cgit v1.2.3