From 44861dcbfeee041223c4aac1ee075e92fa4daa01 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sun, 18 Sep 2016 12:42:15 +0200 Subject: update --- eigen/test/sparse_permutations.cpp | 187 +++++++++++++++++++++++++++++++++++++ 1 file changed, 187 insertions(+) create mode 100644 eigen/test/sparse_permutations.cpp (limited to 'eigen/test/sparse_permutations.cpp') diff --git a/eigen/test/sparse_permutations.cpp b/eigen/test/sparse_permutations.cpp new file mode 100644 index 0000000..e4ce1d6 --- /dev/null +++ b/eigen/test/sparse_permutations.cpp @@ -0,0 +1,187 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 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 "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 typename SparseMatrixType::Index Index; + typedef SparseMatrix OtherSparseMatrixType; + 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); + OtherSparseMatrixType 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*p; + res_d = mat_d*p; + VERIFY(res.isApprox(res_d) && "mat*p"); + + res = p*mat; + res_d = p*mat_d; + VERIFY(res.isApprox(res_d) && "p*mat"); + + res = mat*p.inverse(); + res_d = mat*p.inverse(); + VERIFY(res.isApprox(res_d) && "mat*inv(p)"); + + res = p.inverse()*mat; + res_d = p.inverse()*mat_d; + VERIFY(res.isApprox(res_d) && "inv(p)*mat"); + + res = mat.twistedBy(p); + res_d = (p * mat_d) * p.inverse(); + VERIFY(res.isApprox(res_d) && "p*mat*inv(p)"); + + + 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 = mat.template selfadjointView(); + res_d = up_sym_d; + VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full"); + + res = mat.template selfadjointView(); + res_d = lo_sym_d; + VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full"); + + res = up.template selfadjointView(); + res_d = up_sym_d; + VERIFY(res.isApprox(res_d) && "upper selfadjoint to full"); + + res = lo.template selfadjointView(); + res_d = lo_sym_d; + VERIFY(res.isApprox(res_d) && "lower selfadjoint full"); + + + res.template selfadjointView() = mat.template selfadjointView(); + res_d = up_sym_d.template triangularView(); + VERIFY(res.isApprox(res_d) && "full selfadjoint upper to upper"); + + res.template selfadjointView() = mat.template selfadjointView(); + res_d = up_sym_d.template triangularView(); + VERIFY(res.isApprox(res_d) && "full selfadjoint upper to lower"); + + res.template selfadjointView() = mat.template selfadjointView(); + res_d = lo_sym_d.template triangularView(); + VERIFY(res.isApprox(res_d) && "full selfadjoint lower to upper"); + + res.template selfadjointView() = mat.template selfadjointView(); + res_d = lo_sym_d.template triangularView(); + VERIFY(res.isApprox(res_d) && "full selfadjoint lower 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 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"); +} + +template void sparse_permutations_all(int size) +{ + CALL_SUBTEST(( sparse_permutations(SparseMatrix(size,size)) )); + CALL_SUBTEST(( sparse_permutations(SparseMatrix(size,size)) )); + CALL_SUBTEST(( sparse_permutations(SparseMatrix(size,size)) )); + CALL_SUBTEST(( sparse_permutations(SparseMatrix(size,size)) )); +} + +void test_sparse_permutations() +{ + for(int i = 0; i < g_repeat; i++) { + int s = Eigen::internal::random(1,50); + CALL_SUBTEST_1(( sparse_permutations_all(s) )); + CALL_SUBTEST_2(( sparse_permutations_all >(s) )); + } +} -- cgit v1.2.3