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 +// +// 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 . + +#ifndef EIGEN_SPARSE_PERMUTATION_H +#define EIGEN_SPARSE_PERMUTATION_H + +// This file implements sparse * permutation products + +namespace internal { + +template +struct traits > +{ + typedef typename remove_all::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, + SparseMatrix >::type ReturnType; +}; + +template +struct permut_sparsematrix_product_retval + : public ReturnByValue > +{ + typedef typename remove_all::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 inline void evalTo(Dest& dst) const + { + if(MoveOuter) + { + SparseMatrix tmp(m_matrix.rows(), m_matrix.cols()); + VectorXi sizes(m_matrix.outerSize()); + for(Index j=0; j tmp(m_matrix.rows(), m_matrix.cols()); + VectorXi sizes(tmp.outerSize()); + sizes.setZero(); + PermutationMatrix perm; + if((Side==OnTheLeft) ^ Transposed) + perm = m_permutation; + else + perm = m_permutation.transpose(); + + for(Index j=0; j +inline const internal::permut_sparsematrix_product_retval, SparseDerived, OnTheRight, false> +operator*(const SparseMatrixBase& matrix, const PermutationBase& perm) +{ + return internal::permut_sparsematrix_product_retval, SparseDerived, OnTheRight, false>(perm, matrix.derived()); +} + +/** \returns the matrix with the permutation applied to the rows + */ +template +inline const internal::permut_sparsematrix_product_retval, SparseDerived, OnTheLeft, false> +operator*( const PermutationBase& perm, const SparseMatrixBase& matrix) +{ + return internal::permut_sparsematrix_product_retval, SparseDerived, OnTheLeft, false>(perm, matrix.derived()); +} + + + +/** \returns the matrix with the inverse permutation applied to the columns. + */ +template +inline const internal::permut_sparsematrix_product_retval, SparseDerived, OnTheRight, true> +operator*(const SparseMatrixBase& matrix, const Transpose >& tperm) +{ + return internal::permut_sparsematrix_product_retval, SparseDerived, OnTheRight, true>(tperm.nestedPermutation(), matrix.derived()); +} + +/** \returns the matrix with the inverse permutation applied to the rows. + */ +template +inline const internal::permut_sparsematrix_product_retval, SparseDerived, OnTheLeft, true> +operator*(const Transpose >& tperm, const SparseMatrixBase& matrix) +{ + return internal::permut_sparsematrix_product_retval, 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 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().twistedBy(p_null); res_d = up_sym_d; @@ -76,6 +96,41 @@ template void sparse_permutations(c 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);