Eigen/src/SparseCore/SparsePermutation.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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/.

#ifndef EIGEN_SPARSE_PERMUTATION_H
#define EIGEN_SPARSE_PERMUTATION_H

// This file implements sparse * permutation products

#include "./InternalHeaderCheck.h"

namespace Eigen { 

namespace internal {

template<typename ExpressionType, typename PlainObjectType, bool NeedEval = !is_same<ExpressionType, PlainObjectType>::value>
struct XprHelper
{
    XprHelper(const ExpressionType& xpr) : m_xpr(xpr) {}
    inline const PlainObjectType& xpr() const { return m_xpr; }
    // this is a new PlainObjectType initialized by xpr
    const PlainObjectType m_xpr;
};
template<typename ExpressionType, typename PlainObjectType>
struct XprHelper<ExpressionType, PlainObjectType, false>
{
    XprHelper(const ExpressionType& xpr) : m_xpr(xpr) {}
    inline const PlainObjectType& xpr() const { return m_xpr; }
    // this is a reference to xpr
    const PlainObjectType& m_xpr;
};

template<typename PermDerived, bool NeedInverseEval>
struct PermHelper
{
    using IndicesType = typename PermDerived::IndicesType;
    using PermutationIndex = typename IndicesType::Scalar;
    using type = PermutationMatrix<IndicesType::SizeAtCompileTime, IndicesType::MaxSizeAtCompileTime, PermutationIndex>;
    PermHelper(const PermDerived& perm) : m_perm(perm.inverse()) {}
    inline const type& perm() const { return m_perm; }
    // this is a new PermutationMatrix initialized by perm.inverse()
    const type m_perm;
};
template<typename PermDerived>
struct PermHelper<PermDerived, false>
{
    using type = PermDerived;
    PermHelper(const PermDerived& perm) : m_perm(perm) {}
    inline const type& perm() const { return m_perm; }
    // this is a reference to perm
    const type& m_perm;
};

template<typename ExpressionType, int Side, bool Transposed>
struct permutation_matrix_product<ExpressionType, Side, Transposed, SparseShape>
{
    using MatrixType = typename nested_eval<ExpressionType, 1>::type;
    using MatrixTypeCleaned = remove_all_t<MatrixType>;

    using Scalar = typename MatrixTypeCleaned::Scalar;
    using StorageIndex = typename MatrixTypeCleaned::StorageIndex;

    // the actual "return type" is `Dest`. this is a temporary type
    using ReturnType = SparseMatrix<Scalar, MatrixTypeCleaned::IsRowMajor ? RowMajor : ColMajor, StorageIndex>;
    using TmpHelper = XprHelper<ExpressionType, ReturnType>;

    static constexpr bool NeedOuterPermutation = ExpressionType::IsRowMajor ? Side == OnTheLeft : Side == OnTheRight;
    static constexpr bool NeedInversePermutation = Transposed ? Side == OnTheLeft : Side == OnTheRight;

    template <typename Dest, typename PermutationType>
    static inline void permute_outer(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) {

        // if ExpressionType is not ReturnType, evaluate `xpr` (allocation)
        // otherwise, just reference `xpr`
        // TODO: handle trivial expressions such as CwiseBinaryOp without temporary
        const TmpHelper tmpHelper(xpr);
        const ReturnType& tmp = tmpHelper.xpr();

        ReturnType result(tmp.rows(), tmp.cols());

        for (Index j = 0; j < tmp.outerSize(); j++) {
          Index jp = perm.indices().coeff(j);
          Index jsrc = NeedInversePermutation ? jp : j;
          Index jdst = NeedInversePermutation ? j : jp;
          Index begin = tmp.outerIndexPtr()[jsrc];
          Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[jsrc + 1] : begin + tmp.innerNonZeroPtr()[jsrc];
          result.outerIndexPtr()[jdst + 1] += end - begin;
        }

        std::partial_sum(result.outerIndexPtr(), result.outerIndexPtr() + result.outerSize() + 1,
                         result.outerIndexPtr());
        result.resizeNonZeros(result.nonZeros());

        for (Index j = 0; j < tmp.outerSize(); j++) {
          Index jp = perm.indices().coeff(j);
          Index jsrc = NeedInversePermutation ? jp : j;
          Index jdst = NeedInversePermutation ? j : jp;
          Index begin = tmp.outerIndexPtr()[jsrc];
          Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[jsrc + 1] : begin + tmp.innerNonZeroPtr()[jsrc];
          Index target = result.outerIndexPtr()[jdst];
          smart_copy(tmp.innerIndexPtr() + begin, tmp.innerIndexPtr() + end, result.innerIndexPtr() + target);
          smart_copy(tmp.valuePtr() + begin, tmp.valuePtr() + end, result.valuePtr() + target);
        }
        dst = std::move(result);
    }

    template <typename Dest, typename PermutationType>
    static inline void permute_inner(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) {
        using InnerPermHelper = PermHelper<PermutationType, NeedInversePermutation>;
        using InnerPermType = typename InnerPermHelper::type;

        // if ExpressionType is not ReturnType, evaluate `xpr` (allocation)
        // otherwise, just reference `xpr`
        // TODO: handle trivial expressions such as CwiseBinaryOp without temporary
        const TmpHelper tmpHelper(xpr);
        const ReturnType& tmp = tmpHelper.xpr();

        // if inverse permutation of inner indices is requested, calculate perm.inverse() (allocation)
        // otherwise, just reference `perm`
        const InnerPermHelper permHelper(perm);
        const InnerPermType& innerPerm = permHelper.perm();

        ReturnType result(tmp.rows(), tmp.cols());

        for (Index j = 0; j < tmp.outerSize(); j++) {
            Index begin = tmp.outerIndexPtr()[j];
            Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[j + 1] : begin + tmp.innerNonZeroPtr()[j];
            result.outerIndexPtr()[j + 1] += end - begin;
        }

        std::partial_sum(result.outerIndexPtr(), result.outerIndexPtr() + result.outerSize() + 1, result.outerIndexPtr());
        result.resizeNonZeros(result.nonZeros());

        for (Index j = 0; j < tmp.outerSize(); j++) {
            Index begin = tmp.outerIndexPtr()[j];
            Index end = tmp.isCompressed() ? tmp.outerIndexPtr()[j + 1] : begin + tmp.innerNonZeroPtr()[j];
            Index target = result.outerIndexPtr()[j];
            std::transform(tmp.innerIndexPtr() + begin, tmp.innerIndexPtr() + end, result.innerIndexPtr() + target,
                           [&innerPerm](StorageIndex i) { return innerPerm.indices().coeff(i); });
            smart_copy(tmp.valuePtr() + begin, tmp.valuePtr() + end, result.valuePtr() + target);
        }
        // the inner indices were permuted, and must be sorted
        result.sortInnerIndices();
        dst = std::move(result);
    }

    template <typename Dest, typename PermutationType, bool DoOuter = NeedOuterPermutation, std::enable_if_t<DoOuter, int> = 0>
    static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { permute_outer(dst, perm, xpr); }

    template <typename Dest, typename PermutationType, bool DoOuter = NeedOuterPermutation, std::enable_if_t<!DoOuter, int> = 0>
    static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) { permute_inner(dst, perm, xpr); }
};

}

namespace internal {

template <int ProductTag> struct product_promote_storage_type<Sparse,             PermutationStorage, ProductTag> { typedef Sparse ret; };
template <int ProductTag> struct product_promote_storage_type<PermutationStorage, Sparse,             ProductTag> { typedef Sparse ret; };

// TODO, the following two overloads are only needed to define the right temporary type through 
// typename traits<permutation_sparse_matrix_product<Rhs,Lhs,OnTheRight,false> >::ReturnType
// whereas it should be correctly handled by traits<Product<> >::PlainObject

template<typename Lhs, typename Rhs, int ProductTag>
struct product_evaluator<Product<Lhs, Rhs, AliasFreeProduct>, ProductTag, PermutationShape, SparseShape>
  : public evaluator<typename permutation_matrix_product<Rhs,OnTheLeft,false,SparseShape>::ReturnType>
{
  typedef Product<Lhs, Rhs, AliasFreeProduct> XprType;
  typedef typename permutation_matrix_product<Rhs,OnTheLeft,false,SparseShape>::ReturnType PlainObject;
  typedef evaluator<PlainObject> Base;

  enum {
    Flags = Base::Flags | EvalBeforeNestingBit
  };

  explicit product_evaluator(const XprType& xpr)
    : m_result(xpr.rows(), xpr.cols())
  {
    internal::construct_at<Base>(this, m_result);
    generic_product_impl<Lhs, Rhs, PermutationShape, SparseShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
  }

protected:
  PlainObject m_result;
};

template<typename Lhs, typename Rhs, int ProductTag>
struct product_evaluator<Product<Lhs, Rhs, AliasFreeProduct>, ProductTag, SparseShape, PermutationShape >
  : public evaluator<typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType>
{
  typedef Product<Lhs, Rhs, AliasFreeProduct> XprType;
  typedef typename permutation_matrix_product<Lhs,OnTheRight,false,SparseShape>::ReturnType PlainObject;
  typedef evaluator<PlainObject> Base;

  enum {
    Flags = Base::Flags | EvalBeforeNestingBit
  };

  explicit product_evaluator(const XprType& xpr)
    : m_result(xpr.rows(), xpr.cols())
  {
    ::new (static_cast<Base*>(this)) Base(m_result);
    generic_product_impl<Lhs, Rhs, SparseShape, PermutationShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
  }

protected:
  PlainObject m_result;
};

} // end namespace internal

/** \returns the matrix with the permutation applied to the columns
 */
template <typename SparseDerived, typename PermDerived>
inline const Product<SparseDerived, PermDerived, AliasFreeProduct> operator*(
    const SparseMatrixBase<SparseDerived>& matrix, const PermutationBase<PermDerived>& perm) {
  return Product<SparseDerived, PermDerived, AliasFreeProduct>(matrix.derived(), perm.derived());
}

/** \returns the matrix with the permutation applied to the rows
 */
template <typename SparseDerived, typename PermDerived>
inline const Product<PermDerived, SparseDerived, AliasFreeProduct> operator*(
    const PermutationBase<PermDerived>& perm, const SparseMatrixBase<SparseDerived>& matrix) {
  return Product<PermDerived, SparseDerived, AliasFreeProduct>(perm.derived(), matrix.derived());
}

/** \returns the matrix with the inverse permutation applied to the columns.
 */
template <typename SparseDerived, typename PermutationType>
inline const Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct> operator*(
    const SparseMatrixBase<SparseDerived>& matrix, const InverseImpl<PermutationType, PermutationStorage>& tperm) {
  return Product<SparseDerived, Inverse<PermutationType>, AliasFreeProduct>(matrix.derived(), tperm.derived());
}

/** \returns the matrix with the inverse permutation applied to the rows.
 */
template <typename SparseDerived, typename PermutationType>
inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct> operator*(
    const InverseImpl<PermutationType, PermutationStorage>& tperm, const SparseMatrixBase<SparseDerived>& matrix) {
  return Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct>(tperm.derived(), matrix.derived());
}

} // end namespace Eigen

#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H