SparsePermutation.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd 
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_SPARSE_PERMUTATION_H
11#define EIGEN_SPARSE_PERMUTATION_H
12
13// This file implements sparse * permutation products
14
15namespace Eigen {
16
17namespace internal {
18
19template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
20struct traits<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
21{
22  typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
23  typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
24  typedef typename MatrixTypeNestedCleaned::Index Index;
25  enum {
26    SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
27    MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
28  };
29
30  typedef typename internal::conditional<MoveOuter,
31        SparseMatrix<Scalar,SrcStorageOrder,Index>,
32        SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> >::type ReturnType;
33};
34
35template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
36struct permut_sparsematrix_product_retval
37 : public ReturnByValue<permut_sparsematrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
38{
39    typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
40    typedef typename MatrixTypeNestedCleaned::Scalar Scalar;
41    typedef typename MatrixTypeNestedCleaned::Index Index;
42
43    enum {
44      SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor,
45      MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight
46    };
47
48    permut_sparsematrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
49      : m_permutation(perm), m_matrix(matrix)
50    {}
51
52    inline int rows() const { return m_matrix.rows(); }
53    inline int cols() const { return m_matrix.cols(); }
54
55    template<typename Dest> inline void evalTo(Dest& dst) const
56    {
57      if(MoveOuter)
58      {
59        SparseMatrix<Scalar,SrcStorageOrder,Index> tmp(m_matrix.rows(), m_matrix.cols());
60        VectorXi sizes(m_matrix.outerSize());
61        for(Index j=0; j<m_matrix.outerSize(); ++j)
62        {
63          Index jp = m_permutation.indices().coeff(j);
64          sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).size();
65        }
66        tmp.reserve(sizes);
67        for(Index j=0; j<m_matrix.outerSize(); ++j)
68        {
69          Index jp = m_permutation.indices().coeff(j);
70          Index jsrc = ((Side==OnTheRight) ^ Transposed) ? jp : j;
71          Index jdst = ((Side==OnTheLeft) ^ Transposed) ? jp : j;
72          for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,jsrc); it; ++it)
73            tmp.insertByOuterInner(jdst,it.index()) = it.value();
74        }
75        dst = tmp;
76      }
77      else
78      {
79        SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> tmp(m_matrix.rows(), m_matrix.cols());
80        VectorXi sizes(tmp.outerSize());
81        sizes.setZero();
82        PermutationMatrix<Dynamic,Dynamic,Index> perm;
83        if((Side==OnTheLeft) ^ Transposed)
84          perm = m_permutation;
85        else
86          perm = m_permutation.transpose();
87
88        for(Index j=0; j<m_matrix.outerSize(); ++j)
89          for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it)
90            sizes[perm.indices().coeff(it.index())]++;
91        tmp.reserve(sizes);
92        for(Index j=0; j<m_matrix.outerSize(); ++j)
93          for(typename MatrixTypeNestedCleaned::InnerIterator it(m_matrix,j); it; ++it)
94            tmp.insertByOuterInner(perm.indices().coeff(it.index()),j) = it.value();
95        dst = tmp;
96      }
97    }
98
99  protected:
100    const PermutationType& m_permutation;
101    typename MatrixType::Nested m_matrix;
102};
103
104}
105
106
107
108/** \returns the matrix with the permutation applied to the columns
109  */
110template<typename SparseDerived, typename PermDerived>
111inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false>
112operator*(const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm)
113{
114  return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, false>(perm, matrix.derived());
115}
116
117/** \returns the matrix with the permutation applied to the rows
118  */
119template<typename SparseDerived, typename PermDerived>
120inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false>
121operator*( const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix)
122{
123  return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, false>(perm, matrix.derived());
124}
125
126
127
128/** \returns the matrix with the inverse permutation applied to the columns.
129  */
130template<typename SparseDerived, typename PermDerived>
131inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true>
132operator*(const SparseMatrixBase<SparseDerived>& matrix, const Transpose<PermutationBase<PermDerived> >& tperm)
133{
134  return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheRight, true>(tperm.nestedPermutation(), matrix.derived());
135}
136
137/** \returns the matrix with the inverse permutation applied to the rows.
138  */
139template<typename SparseDerived, typename PermDerived>
140inline const internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true>
141operator*(const Transpose<PermutationBase<PermDerived> >& tperm, const SparseMatrixBase<SparseDerived>& matrix)
142{
143  return internal::permut_sparsematrix_product_retval<PermutationBase<PermDerived>, SparseDerived, OnTheLeft, true>(tperm.nestedPermutation(), matrix.derived());
144}
145
146} // end namespace Eigen
147
148#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H
149