diff options
| -rw-r--r-- | Eigen/SparseCore | 1 | ||||
| -rw-r--r-- | Eigen/src/SparseCore/SparsePermutation.h | 160 | ||||
| -rw-r--r-- | test/sparse_permutations.cpp | 55 |
3 files changed, 216 insertions, 0 deletions
diff --git a/Eigen/SparseCore b/Eigen/SparseCore index 9a7d11124..3d355a90e 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore @@ -44,6 +44,7 @@ struct Sparse {}; #include "src/SparseCore/SparseCwiseUnaryOp.h" #include "src/SparseCore/SparseCwiseBinaryOp.h" #include "src/SparseCore/SparseDot.h" +#include "src/SparseCore/SparsePermutation.h" #include "src/SparseCore/SparseAssign.h" #include "src/SparseCore/SparseRedux.h" #include "src/SparseCore/SparseFuzzy.h" diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h new file mode 100644 index 000000000..1fe2f404a --- /dev/null +++ b/Eigen/src/SparseCore/SparsePermutation.h @@ -0,0 +1,160 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// 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 <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_PERMUTATION_H +#define EIGEN_SPARSE_PERMUTATION_H + +// This file implements sparse * permutation products + +namespace internal { + +template<typename PermutationType, typename MatrixType, int Side, bool Transposed> +struct traits<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> > +{ + typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned; + typedef typename MatrixTypeNestedCleaned::Scalar Scalar; + typedef typename MatrixTypeNestedCleaned::Index Index; + enum { + SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor, + MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight + }; + + typedef typename internal::conditional<MoveOuter, + SparseMatrix<Scalar,SrcStorageOrder,Index>, + SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> >::type ReturnType; +}; + +template<typename PermutationType, typename MatrixType, int Side, bool Transposed> +struct permut_sparsematrix_product_retval + : public ReturnByValue<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> > +{ + typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned; + typedef typename MatrixTypeNestedCleaned::Scalar Scalar; + typedef typename MatrixTypeNestedCleaned::Index Index; + + enum { + SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor, + MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight + }; + + permut_sparsematrix_product_retval(const PermutationType& perm, const MatrixType& matrix) + : m_permutation(perm), m_matrix(matrix) + {} + + inline int rows() const { return m_matrix.rows(); } + inline int cols() const { return m_matrix.cols(); } + + template<typename Dest> inline void evalTo(Dest& dst) const + { + if(MoveOuter) + { + SparseMatrix<Scalar,SrcStorageOrder,Index> tmp(m_matrix.rows(), m_matrix.cols()); + VectorXi sizes(m_matrix.outerSize()); + for(Index j=0; j<m_matrix.outerSize(); ++j) + { + Index jp = m_permutation.indices().coeff(j); + sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).size(); + } + tmp.reserve(sizes); + for(Index j=0; j<m_matrix.outerSize(); ++j) + { + Index jp = m_permutation.indices().coeff(j); + Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j; + Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j; + for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,jsrc); it; ++it) + tmp.insertByOuterInner(jdst,it.index()) = it.value(); + } + dst = tmp; + } + else + { + SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> tmp(m_matrix.rows(), m_matrix.cols()); + VectorXi sizes(tmp.outerSize()); + sizes.setZero(); + PermutationMatrix<Dynamic,Dynamic,Index> perm; + if((Side==OnTheLeft) ^ Transposed) + perm = m_permutation; + else + perm = m_permutation.transpose(); + + for(Index j=0; j<m_matrix.outerSize(); ++j) + for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it) + sizes[perm.indices().coeff(it.index())]++; + tmp.reserve(sizes); + for(Index j=0; j<m_matrix.outerSize(); ++j) + for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it) + tmp.insertByOuterInner(perm.indices().coeff(it.index()),j) = it.value(); + dst = tmp; + } + } + + protected: + const PermutationType& m_permutation; + typename MatrixType::Nested m_matrix; +}; + +} + + + +/** \returns the matrix with the permutation applied to the columns + */ +template<typename SparseDerived, typename PermDerived> +inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false> +operator*(const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm) +{ + return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false>(perm, matrix.derived()); +} + +/** \returns the matrix with the permutation applied to the rows + */ +template<typename SparseDerived, typename PermDerived> +inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false> +operator*( const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix) +{ + return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false>(perm, matrix.derived()); +} + + + +/** \returns the matrix with the inverse permutation applied to the columns. + */ +template<typename SparseDerived, typename PermDerived> +inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true> +operator*(const SparseMatrixBase<SparseDerived>& matrix, const Transpose<PermutationBase<PermDerived> >& tperm) +{ + return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true>(tperm.nestedPermutation(), matrix.derived()); +} + +/** \returns the matrix with the inverse permutation applied to the rows. + */ +template<typename SparseDerived, typename PermDerived> +inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true> +operator*(const Transpose<PermutationBase<PermDerived> >& tperm, const SparseMatrixBase<SparseDerived>& matrix) +{ + return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true>(tperm.nestedPermutation(), matrix.derived()); +} + + +#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H diff --git a/test/sparse_permutations.cpp b/test/sparse_permutations.cpp index 522e78f67..4d22e358e 100644 --- a/test/sparse_permutations.cpp +++ b/test/sparse_permutations.cpp @@ -59,6 +59,26 @@ template<int OtherStorage, typename SparseMatrixType> void sparse_permutations(c 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<Upper>().twistedBy(p_null); res_d = up_sym_d; @@ -76,6 +96,41 @@ template<int OtherStorage, typename SparseMatrixType> void sparse_permutations(c res = lo.template selfadjointView<Lower>().twistedBy(p_null); res_d = lo_sym_d; VERIFY(res.isApprox(res_d) && "lower selfadjoint full"); + + + res = mat.template selfadjointView<Upper>(); + res_d = up_sym_d; + VERIFY(res.isApprox(res_d) && "full selfadjoint upper to full"); + + res = mat.template selfadjointView<Lower>(); + res_d = lo_sym_d; + VERIFY(res.isApprox(res_d) && "full selfadjoint lower to full"); + + res = up.template selfadjointView<Upper>(); + res_d = up_sym_d; + VERIFY(res.isApprox(res_d) && "upper selfadjoint to full"); + + res = lo.template selfadjointView<Lower>(); + res_d = lo_sym_d; + VERIFY(res.isApprox(res_d) && "lower selfadjoint full"); + + + res.template selfadjointView<Upper>() = mat.template selfadjointView<Upper>(); + res_d = up_sym_d.template triangularView<Upper>(); + VERIFY(res.isApprox(res_d) && "full selfadjoint upper to upper"); + + res.template selfadjointView<Lower>() = mat.template selfadjointView<Upper>(); + res_d = up_sym_d.template triangularView<Lower>(); + VERIFY(res.isApprox(res_d) && "full selfadjoint upper to lower"); + + res.template selfadjointView<Upper>() = mat.template selfadjointView<Lower>(); + res_d = lo_sym_d.template triangularView<Upper>(); + VERIFY(res.isApprox(res_d) && "full selfadjoint lower to upper"); + + res.template selfadjointView<Lower>() = mat.template selfadjointView<Lower>(); + res_d = lo_sym_d.template triangularView<Lower>(); + VERIFY(res.isApprox(res_d) && "full selfadjoint lower to lower"); + res.template selfadjointView<Upper>() = mat.template selfadjointView<Upper>().twistedBy(p); |