ScaledPermutationMatrix.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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/.
// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_SCALEDPERMUTATIONMATRIX_H
#define EIGEN_SCALEDPERMUTATIONMATRIX_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
namespace internal {
 
template <typename Scalar_, int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_>
struct traits<ScaledPermutationMatrix<Scalar_, SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>> {
  using Scalar = Scalar_;
  using StorageIndex = StorageIndex_;
  // Dense, so that products and sums with dense operands promote with the existing rules; the O(n) kernels are
  // selected by the shape (ScaledPermutationShape), as for TriangularView.
  using StorageKind = Dense;
  using XprKind = MatrixXpr;
  using PermutationType = PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>;
  using ScalesType = Matrix<Scalar_, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>;
  static constexpr int RowsAtCompileTime = SizeAtCompileTime;
  static constexpr int ColsAtCompileTime = SizeAtCompileTime;
  static constexpr int MaxRowsAtCompileTime = MaxSizeAtCompileTime;
  static constexpr int MaxColsAtCompileTime = MaxSizeAtCompileTime;
  static constexpr unsigned int Flags = NestByRefBit;
};

template <typename Scalar, typename Lhs, typename Rhs, typename StorageIndex>
struct scaled_permutation_result {
  static constexpr int Size = size_prefer_fixed(traits<Lhs>::RowsAtCompileTime, traits<Rhs>::RowsAtCompileTime);
  static constexpr int MaxSize =
      min_size_prefer_fixed(traits<Lhs>::MaxRowsAtCompileTime, traits<Rhs>::MaxRowsAtCompileTime);
  using type = ScaledPermutationMatrix<Scalar, Size, MaxSize, StorageIndex>;
};
 
template <typename Lhs, typename Rhs, typename StorageIndex>
struct scaled_permutation_product_result {
  using LhsScalar = typename traits<Lhs>::Scalar;
  using RhsScalar = typename traits<Rhs>::Scalar;
  using Scalar =
      typename ScalarBinaryOpTraits<LhsScalar, RhsScalar, scalar_product_op<LhsScalar, RhsScalar>>::ReturnType;
  using type = typename scaled_permutation_result<Scalar, Lhs, Rhs, StorageIndex>::type;
};

template <typename IndicesType, typename Result>
EIGEN_DEVICE_FUNC void invert_permutation_indices(const IndicesType& indices, Result& result) {
  using StorageIndex = typename Result::Scalar;
  for (Index k = 0; k < indices.size(); ++k) result.coeffRef(indices.coeff(k)) = StorageIndex(k);
}
 
}  // namespace internal

template <typename Derived>
class ScaledPermutationBase : public EigenBase<Derived> {
  using Traits = internal::traits<Derived>;
  using Base = EigenBase<Derived>;
 
 public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
  using Scalar = typename Traits::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using StorageIndex = typename Traits::StorageIndex;
  using StorageKind = typename Traits::StorageKind;
  using PermutationType = typename Traits::PermutationType;
  using ScalesType = typename Traits::ScalesType;
  using IndicesType = typename PermutationType::IndicesType;
  static constexpr int RowsAtCompileTime = Traits::RowsAtCompileTime;
  static constexpr int ColsAtCompileTime = Traits::ColsAtCompileTime;
  static constexpr int MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime;
  static constexpr int MaxColsAtCompileTime = Traits::MaxColsAtCompileTime;
  static constexpr int SizeAtCompileTime = internal::size_at_compile_time(RowsAtCompileTime, ColsAtCompileTime);
  static constexpr int MaxSizeAtCompileTime =
      internal::size_at_compile_time(MaxRowsAtCompileTime, MaxColsAtCompileTime);
  static constexpr bool IsVectorAtCompileTime = false;
  static constexpr unsigned int Flags = Traits::Flags;
  using DenseMatrixType =
      Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime>;
  using PlainObject = ScaledPermutationMatrix<Scalar, RowsAtCompileTime, MaxRowsAtCompileTime, StorageIndex>;
  using Nested = const Derived&;
  using Base::derived;
#endif

  EIGEN_DEVICE_FUNC constexpr Index rows() const { return indices().size(); }
  EIGEN_DEVICE_FUNC constexpr Index cols() const { return indices().size(); }

  EIGEN_DEVICE_FUNC const PermutationType& permutation() const { return derived().permutation(); }
  EIGEN_DEVICE_FUNC constexpr const IndicesType& indices() const { return derived().indices(); }
  EIGEN_DEVICE_FUNC const ScalesType& scales() const { return derived().scales(); }

  EIGEN_DEVICE_FUNC Scalar coeff(Index row, Index col) const {
    eigen_assert(row >= 0 && col >= 0 && row < rows() && col < cols());
    return Index(indices().coeff(col)) == row ? scales().coeff(col) : Scalar(0);
  }

  EIGEN_DEVICE_FUNC DenseMatrixType toDenseMatrix() const { return derived(); }
 
#ifndef EIGEN_PARSED_BY_DOXYGEN
  template <typename DenseDerived>
  EIGEN_DEVICE_FUNC void evalTo(MatrixBase<DenseDerived>& other) const {
    other.setZero();
    for (Index k = 0; k < rows(); ++k) other.coeffRef(indices().coeff(k), k) = scales().coeff(k);
  }
#endif

  PlainObject inverse() const {
    PlainObject result(rows());
    for (Index k = 0; k < rows(); ++k) {
      const Index image = indices().coeff(k);
      result.indices().coeffRef(image) = StorageIndex(k);
      result.scales().coeffRef(image) = Scalar(1) / scales().coeff(k);
    }
    return result;
  }

  PlainObject transpose() const {
    PlainObject result(rows());
    for (Index k = 0; k < rows(); ++k) {
      const Index image = indices().coeff(k);
      result.indices().coeffRef(image) = StorageIndex(k);
      result.scales().coeffRef(image) = scales().coeff(k);
    }
    return result;
  }

  PlainObject adjoint() const {
    PlainObject result = transpose();
    result.scales() = result.scales().conjugate();
    return result;
  }

  Scalar determinant() const { return Scalar(RealScalar(permutation().determinant())) * scales().prod(); }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC const Product<Derived, OtherDerived, DefaultProduct> operator*(
      const MatrixBase<OtherDerived>& other) const {
    return Product<Derived, OtherDerived, DefaultProduct>(derived(), other.derived());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC friend const Product<OtherDerived, Derived, DefaultProduct> operator*(
      const MatrixBase<OtherDerived>& other, const ScaledPermutationBase& scaled) {
    return Product<OtherDerived, Derived, DefaultProduct>(other.derived(), scaled.derived());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC const Product<Derived, OtherDerived, DefaultProduct> operator*(
      const TriangularBase<OtherDerived>& other) const {
    return Product<Derived, OtherDerived, DefaultProduct>(derived(), other.derived());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC friend const Product<OtherDerived, Derived, DefaultProduct> operator*(
      const TriangularBase<OtherDerived>& other, const ScaledPermutationBase& scaled) {
    return Product<OtherDerived, Derived, DefaultProduct>(other.derived(), scaled.derived());
  }
 
  // Closed-form products, all O(n). With S e_c = v_c e_{tau(c)} the rules are:
  //   S D: scale v_c d_c            D S: scale d_{tau(c)} v_c
  //   S P: (tau(sigma(c)), v_{sigma(c)})   P S: (sigma(tau(c)), v_c)
  //   S T: (tau(tau2(c)), v_{tau2(c)} v2_c)

  template <typename OtherDerived>
  typename internal::scaled_permutation_product_result<Derived, OtherDerived, StorageIndex>::type operator*(
      const DiagonalBase<OtherDerived>& other) const {
    using Result = typename internal::scaled_permutation_product_result<Derived, OtherDerived, StorageIndex>::type;
    eigen_assert(rows() == other.rows());
    return Result(indices(), scales().cwiseProduct(other.diagonal()));
  }

  template <typename OtherDerived>
  friend typename internal::scaled_permutation_product_result<OtherDerived, Derived, StorageIndex>::type operator*(
      const DiagonalBase<OtherDerived>& other, const ScaledPermutationBase& scaled) {
    using Result = typename internal::scaled_permutation_product_result<OtherDerived, Derived, StorageIndex>::type;
    eigen_assert(other.rows() == scaled.rows());
    Result result(scaled.rows());
    result.indices() = scaled.indices();
    for (Index k = 0; k < scaled.rows(); ++k)
      result.scales().coeffRef(k) = other.diagonal().coeff(scaled.indices().coeff(k)) * scaled.scales().coeff(k);
    return result;
  }

  template <typename OtherDerived>
  PlainObject operator*(const PermutationBase<OtherDerived>& other) const {
    eigen_assert(rows() == other.rows());
    PlainObject result(rows());
    for (Index k = 0; k < rows(); ++k) {
      const Index source = other.indices().coeff(k);
      result.indices().coeffRef(k) = indices().coeff(source);
      result.scales().coeffRef(k) = scales().coeff(source);
    }
    return result;
  }

  template <typename OtherDerived>
  friend PlainObject operator*(const PermutationBase<OtherDerived>& other, const ScaledPermutationBase& scaled) {
    eigen_assert(other.rows() == scaled.rows());
    PlainObject result(scaled.rows());
    for (Index k = 0; k < scaled.rows(); ++k)
      result.indices().coeffRef(k) = StorageIndex(other.indices().coeff(scaled.indices().coeff(k)));
    result.scales() = scaled.scales();
    return result;
  }

  template <typename OtherDerived>
  PlainObject operator*(const InverseImpl<OtherDerived, PermutationStorage>& other) const {
    eigen_assert(rows() == other.rows());
    IndicesType inverse(rows());
    internal::invert_permutation_indices(other.derived().nestedExpression().indices(), inverse);
    PlainObject result(rows());
    for (Index k = 0; k < rows(); ++k) {
      const Index source = inverse.coeff(k);
      result.indices().coeffRef(k) = indices().coeff(source);
      result.scales().coeffRef(k) = scales().coeff(source);
    }
    return result;
  }

  template <typename OtherDerived>
  friend PlainObject operator*(const InverseImpl<OtherDerived, PermutationStorage>& other,
                               const ScaledPermutationBase& scaled) {
    eigen_assert(other.rows() == scaled.rows());
    IndicesType inverse(scaled.rows());
    internal::invert_permutation_indices(other.derived().nestedExpression().indices(), inverse);
    PlainObject result(scaled.rows());
    for (Index k = 0; k < scaled.rows(); ++k) result.indices().coeffRef(k) = inverse.coeff(scaled.indices().coeff(k));
    result.scales() = scaled.scales();
    return result;
  }

  template <typename OtherDerived>
  typename internal::scaled_permutation_product_result<Derived, OtherDerived, StorageIndex>::type operator*(
      const ScaledPermutationBase<OtherDerived>& other) const {
    using Result = typename internal::scaled_permutation_product_result<Derived, OtherDerived, StorageIndex>::type;
    eigen_assert(rows() == other.rows());
    Result result(rows());
    for (Index k = 0; k < rows(); ++k) {
      const Index source = other.indices().coeff(k);
      result.indices().coeffRef(k) = indices().coeff(source);
      result.scales().coeffRef(k) = scales().coeff(source) * other.scales().coeff(k);
    }
    return result;
  }

  PlainObject operator*(const Scalar& alpha) const { return PlainObject(indices(), scales() * alpha); }

  friend PlainObject operator*(const Scalar& alpha, const ScaledPermutationBase& scaled) {
    return PlainObject(scaled.indices(), alpha * scaled.scales());
  }

  PlainObject operator-() const { return PlainObject(indices(), -scales()); }
 
  // Sums with a dense matrix are lazy: the scaled permutation is read through its index-based evaluator.

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC friend const EIGEN_CWISE_BINARY_RETURN_TYPE(OtherDerived, Derived, internal::scalar_sum_op)
  operator+(const MatrixBase<OtherDerived>& lhs, const ScaledPermutationBase & rhs) {
    return EIGEN_CWISE_BINARY_RETURN_TYPE(OtherDerived, Derived, internal::scalar_sum_op)(lhs.derived(), rhs.derived());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC friend const EIGEN_CWISE_BINARY_RETURN_TYPE(Derived, OtherDerived, internal::scalar_sum_op)
  operator+(const ScaledPermutationBase & lhs, const MatrixBase<OtherDerived>& rhs) {
    return EIGEN_CWISE_BINARY_RETURN_TYPE(Derived, OtherDerived, internal::scalar_sum_op)(lhs.derived(), rhs.derived());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC friend const EIGEN_CWISE_BINARY_RETURN_TYPE(OtherDerived, Derived, internal::scalar_difference_op)
  operator-(const MatrixBase<OtherDerived>& lhs, const ScaledPermutationBase & rhs) {
    return EIGEN_CWISE_BINARY_RETURN_TYPE(OtherDerived, Derived, internal::scalar_difference_op)(lhs.derived(),
                                                                                                 rhs.derived());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC friend const EIGEN_CWISE_BINARY_RETURN_TYPE(Derived, OtherDerived, internal::scalar_difference_op)
  operator-(const ScaledPermutationBase & lhs, const MatrixBase<OtherDerived>& rhs) {
    return EIGEN_CWISE_BINARY_RETURN_TYPE(Derived, OtherDerived, internal::scalar_difference_op)(lhs.derived(),
                                                                                                 rhs.derived());
  }
};

template <typename Scalar_, int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_>
class ScaledPermutationMatrix
    : public ScaledPermutationBase<
          ScaledPermutationMatrix<Scalar_, SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>> {
  using Base = ScaledPermutationBase<ScaledPermutationMatrix>;
  using Traits = internal::traits<ScaledPermutationMatrix>;
 
 public:
#ifndef EIGEN_PARSED_BY_DOXYGEN
  using Scalar = Scalar_;
  using StorageIndex = StorageIndex_;
  using PermutationType = typename Traits::PermutationType;
  using ScalesType = typename Traits::ScalesType;
  using IndicesType = typename PermutationType::IndicesType;
#endif

  EIGEN_DEVICE_FUNC ScaledPermutationMatrix() = default;

  EIGEN_DEVICE_FUNC explicit ScaledPermutationMatrix(Index size) : m_permutation(size), m_scales(size) {}

  template <typename PermutationDerived, typename DiagonalDerived>
  EIGEN_DEVICE_FUNC ScaledPermutationMatrix(const PermutationBase<PermutationDerived>& permutation,
                                            const DiagonalBase<DiagonalDerived>& diagonal)
      : m_permutation(permutation), m_scales(diagonal.diagonal()) {
    eigen_assert(permutation.rows() == diagonal.rows());
  }

  template <typename IndicesDerived, typename ScalesDerived>
  EIGEN_DEVICE_FUNC ScaledPermutationMatrix(const MatrixBase<IndicesDerived>& indices,
                                            const MatrixBase<ScalesDerived>& scales)
      : m_permutation(indices), m_scales(scales) {
    eigen_assert(indices.size() == scales.size());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC explicit ScaledPermutationMatrix(const PermutationBase<OtherDerived>& permutation)
      : m_permutation(permutation), m_scales(ScalesType::Ones(permutation.rows())) {}

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC explicit ScaledPermutationMatrix(const DiagonalBase<OtherDerived>& diagonal)
      : m_permutation(diagonal.rows()), m_scales(diagonal.diagonal()) {
    m_permutation.setIdentity();
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC ScaledPermutationMatrix(const ScaledPermutationBase<OtherDerived>& other)
      : m_permutation(other.permutation()), m_scales(other.scales()) {}

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC ScaledPermutationMatrix& operator=(const ScaledPermutationBase<OtherDerived>& other) {
    m_permutation = other.permutation();
    m_scales = other.scales();
    return *this;
  }

  EIGEN_DEVICE_FUNC const PermutationType& permutation() const { return m_permutation; }
  EIGEN_DEVICE_FUNC PermutationType& permutation() { return m_permutation; }
  EIGEN_DEVICE_FUNC constexpr const IndicesType& indices() const { return m_permutation.indices(); }
  EIGEN_DEVICE_FUNC constexpr IndicesType& indices() { return m_permutation.indices(); }
  EIGEN_DEVICE_FUNC const ScalesType& scales() const { return m_scales; }
  EIGEN_DEVICE_FUNC ScalesType& scales() { return m_scales; }

  EIGEN_DEVICE_FUNC void resize(Index size) {
    m_permutation.resize(size);
    m_scales.resize(size);
  }

  EIGEN_DEVICE_FUNC void setIdentity() {
    m_permutation.setIdentity();
    m_scales.setOnes();
  }

  EIGEN_DEVICE_FUNC void setIdentity(Index size) {
    resize(size);
    setIdentity();
  }
 
 protected:
  PermutationType m_permutation;
  ScalesType m_scales;
};
 
// The products of a permutation (or its inverse) and a diagonal matrix. With P e_c = e_{sigma(c)}:
//   P D:   (sigma(c), d_c)          D P:   (sigma(c), d_{sigma(c)})
//   P^-1 D: (sigma^-1(c), d_c)      D P^-1: (sigma^-1(c), d_{sigma^-1(c)})

template <typename PermutationDerived, typename DiagonalDerived>
typename internal::scaled_permutation_result<typename DiagonalDerived::Scalar, PermutationDerived, DiagonalDerived,
                                             typename PermutationDerived::StorageIndex>::type
operator*(const PermutationBase<PermutationDerived>& permutation, const DiagonalBase<DiagonalDerived>& diagonal) {
  using Result =
      typename internal::scaled_permutation_result<typename DiagonalDerived::Scalar, PermutationDerived,
                                                   DiagonalDerived, typename PermutationDerived::StorageIndex>::type;
  eigen_assert(permutation.rows() == diagonal.rows());
  return Result(permutation.indices(), diagonal.diagonal());
}

template <typename DiagonalDerived, typename PermutationDerived>
typename internal::scaled_permutation_result<typename DiagonalDerived::Scalar, PermutationDerived, DiagonalDerived,
                                             typename PermutationDerived::StorageIndex>::type
operator*(const DiagonalBase<DiagonalDerived>& diagonal, const PermutationBase<PermutationDerived>& permutation) {
  using Result =
      typename internal::scaled_permutation_result<typename DiagonalDerived::Scalar, PermutationDerived,
                                                   DiagonalDerived, typename PermutationDerived::StorageIndex>::type;
  eigen_assert(permutation.rows() == diagonal.rows());
  Result result(permutation.rows());
  result.indices() = permutation.indices();
  for (Index k = 0; k < permutation.rows(); ++k)
    result.scales().coeffRef(k) = diagonal.diagonal().coeff(permutation.indices().coeff(k));
  return result;
}

template <typename PermutationType, typename DiagonalDerived>
typename internal::scaled_permutation_result<typename DiagonalDerived::Scalar, PermutationType, DiagonalDerived,
                                             typename PermutationType::StorageIndex>::type
operator*(const InverseImpl<PermutationType, PermutationStorage>& inverse,
          const DiagonalBase<DiagonalDerived>& diagonal) {
  using Result =
      typename internal::scaled_permutation_result<typename DiagonalDerived::Scalar, PermutationType, DiagonalDerived,
                                                   typename PermutationType::StorageIndex>::type;
  eigen_assert(inverse.rows() == diagonal.rows());
  Result result(inverse.rows());
  internal::invert_permutation_indices(inverse.derived().nestedExpression().indices(), result.indices());
  result.scales() = diagonal.diagonal();
  return result;
}

template <typename DiagonalDerived, typename PermutationType>
typename internal::scaled_permutation_result<typename DiagonalDerived::Scalar, PermutationType, DiagonalDerived,
                                             typename PermutationType::StorageIndex>::type
operator*(const DiagonalBase<DiagonalDerived>& diagonal,
          const InverseImpl<PermutationType, PermutationStorage>& inverse) {
  using Result =
      typename internal::scaled_permutation_result<typename DiagonalDerived::Scalar, PermutationType, DiagonalDerived,
                                                   typename PermutationType::StorageIndex>::type;
  eigen_assert(inverse.rows() == diagonal.rows());
  Result result(inverse.rows());
  internal::invert_permutation_indices(inverse.derived().nestedExpression().indices(), result.indices());
  for (Index k = 0; k < inverse.rows(); ++k)
    result.scales().coeffRef(k) = diagonal.diagonal().coeff(result.indices().coeff(k));
  return result;
}
 
namespace internal {
 
template <typename Scalar_, int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_>
struct evaluator_traits<ScaledPermutationMatrix<Scalar_, SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>> {
  using Kind = IndexBased;
  using Shape = ScaledPermutationShape;
};

template <typename Scalar_, int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_>
struct evaluator<ScaledPermutationMatrix<Scalar_, SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>>
    : evaluator_base<ScaledPermutationMatrix<Scalar_, SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>> {
  using XprType = ScaledPermutationMatrix<Scalar_, SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>;
  using Scalar = Scalar_;
  using CoeffReturnType = Scalar;
 
  static constexpr int CoeffReadCost =
      int(NumTraits<StorageIndex_>::ReadCost) + int(NumTraits<Scalar>::ReadCost) + int(NumTraits<Scalar>::AddCost);
  static constexpr unsigned int Flags = 0;
  static constexpr int Alignment = 0;
 
  EIGEN_DEVICE_FUNC explicit evaluator(const XprType& xpr) : m_indices(xpr.indices()), m_scales(xpr.scales()) {}
 
  EIGEN_DEVICE_FUNC Scalar coeff(Index row, Index col) const {
    return Index(m_indices.coeff(col)) == row ? m_scales.coeff(col) : Scalar(0);
  }
 
  // Linear access is requested only for vector-shaped operands (inner products), i.e. a 1x1 matrix.
  EIGEN_DEVICE_FUNC Scalar coeff(Index index) const {
    eigen_assert(index == 0);
    return m_scales.coeff(index);
  }
 
 protected:
  evaluator<typename XprType::IndicesType> m_indices;
  evaluator<typename XprType::ScalesType> m_scales;
};
 
struct ScaledPermutation2Dense {};
 
template <>
struct AssignmentKind<DenseShape, ScaledPermutationShape> {
  using Kind = ScaledPermutation2Dense;
};
 
// Scaled permutation to dense assignment: zero fill plus one scatter of n coefficients.
template <typename DstXprType, typename SrcXprType, typename Functor>
struct Assignment<DstXprType, SrcXprType, Functor, ScaledPermutation2Dense> {
  static EIGEN_DEVICE_FUNC void run(
      DstXprType& dst, const SrcXprType& src,
      const internal::assign_op<typename DstXprType::Scalar, typename SrcXprType::Scalar>&) {
    if (dst.rows() != src.rows() || dst.cols() != src.cols()) dst.resize(src.rows(), src.cols());
    dst.setZero();
    for (Index k = 0; k < src.rows(); ++k) dst.coeffRef(src.indices().coeff(k), k) = src.scales().coeff(k);
  }
 
  static EIGEN_DEVICE_FUNC void run(
      DstXprType& dst, const SrcXprType& src,
      const internal::add_assign_op<typename DstXprType::Scalar, typename SrcXprType::Scalar>&) {
    eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
    for (Index k = 0; k < src.rows(); ++k) dst.coeffRef(src.indices().coeff(k), k) += src.scales().coeff(k);
  }
 
  static EIGEN_DEVICE_FUNC void run(
      DstXprType& dst, const SrcXprType& src,
      const internal::sub_assign_op<typename DstXprType::Scalar, typename SrcXprType::Scalar>&) {
    eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
    for (Index k = 0; k < src.rows(); ++k) dst.coeffRef(src.indices().coeff(k), k) -= src.scales().coeff(k);
  }
};
 
/***************************************************************************
 * Products with a dense matrix: rows (S * X) or columns (X * S) are scaled and permuted, O(n m).
 * evalTo reuses the dense permutation kernel, which handles X aliasing dst, and scales in place afterwards;
 * the accumulating form works on the nested (evaluated once) operand.
 * scaleAndAddTo evaluates dst += alpha * (lhs * rhs) with alpha a left factor of the product. Scalar
 * multiplication need not commute, so both kernels keep that operand order: S * X regroups it as
 * (alpha * v_c) * X(c, j) by associativity; X * S keeps alpha * (X(i, tau(c)) * v_c), as folding alpha into
 * v_c would move it past the dense coefficient.
 ***************************************************************************/
 
template <typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Rhs, ScaledPermutationShape, DenseShape, ProductTag>
    : generic_product_impl_base<Lhs, Rhs,
                                generic_product_impl<Lhs, Rhs, ScaledPermutationShape, DenseShape, ProductTag>> {
  using Scalar = typename Product<Lhs, Rhs>::Scalar;
 
  template <typename Dest>
  static EIGEN_DEVICE_FUNC void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) {
    permutation_matrix_product<Rhs, OnTheLeft, false, DenseShape>::run(dst, lhs.permutation(), rhs);
    for (Index k = 0; k < lhs.rows(); ++k) {
      auto row = dst.row(lhs.indices().coeff(k));
      row = lhs.scales().coeff(k) * row;
    }
  }
 
  template <typename Dest>
  static EIGEN_DEVICE_FUNC void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) {
    typename nested_eval<Rhs, 1>::type rhsNested(rhs);
    for (Index k = 0; k < lhs.rows(); ++k)
      dst.row(lhs.indices().coeff(k)) += (alpha * lhs.scales().coeff(k)) * rhsNested.row(k);
  }
};
 
template <typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Rhs, DenseShape, ScaledPermutationShape, ProductTag>
    : generic_product_impl_base<Lhs, Rhs,
                                generic_product_impl<Lhs, Rhs, DenseShape, ScaledPermutationShape, ProductTag>> {
  using Scalar = typename Product<Lhs, Rhs>::Scalar;
 
  template <typename Dest>
  static EIGEN_DEVICE_FUNC void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) {
    permutation_matrix_product<Lhs, OnTheRight, false, DenseShape>::run(dst, rhs.permutation(), lhs);
    for (Index k = 0; k < rhs.rows(); ++k) dst.col(k) *= rhs.scales().coeff(k);
  }
 
  template <typename Dest>
  static EIGEN_DEVICE_FUNC void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) {
    typename nested_eval<Lhs, 1>::type lhsNested(lhs);
    for (Index k = 0; k < rhs.rows(); ++k)
      dst.col(k) += alpha * (lhsNested.col(rhs.indices().coeff(k)) * rhs.scales().coeff(k));
  }
};
 
// A triangular or self-adjoint operand is evaluated into a plain matrix first, as toDenseMatrix() does.
template <typename Lhs, typename Rhs, int ProductTag, bool ViewOnTheRight>
struct scaled_permutation_view_product_impl;
 
template <typename Lhs, typename Rhs, int ProductTag>
struct scaled_permutation_view_product_impl<Lhs, Rhs, ProductTag, true>
    : generic_product_impl_base<Lhs, Rhs, scaled_permutation_view_product_impl<Lhs, Rhs, ProductTag, true>> {
  using Scalar = typename Product<Lhs, Rhs>::Scalar;
  using DenseType = typename Rhs::DenseMatrixType;
  using DenseImpl = generic_product_impl<Lhs, DenseType, ScaledPermutationShape, DenseShape, ProductTag>;
 
  template <typename Dest>
  static EIGEN_DEVICE_FUNC void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) {
    const DenseType dense(rhs);
    DenseImpl::evalTo(dst, lhs, dense);
  }
  template <typename Dest>
  static EIGEN_DEVICE_FUNC void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) {
    const DenseType dense(rhs);
    DenseImpl::scaleAndAddTo(dst, lhs, dense, alpha);
  }
};
 
template <typename Lhs, typename Rhs, int ProductTag>
struct scaled_permutation_view_product_impl<Lhs, Rhs, ProductTag, false>
    : generic_product_impl_base<Lhs, Rhs, scaled_permutation_view_product_impl<Lhs, Rhs, ProductTag, false>> {
  using Scalar = typename Product<Lhs, Rhs>::Scalar;
  using DenseType = typename Lhs::DenseMatrixType;
  using DenseImpl = generic_product_impl<DenseType, Rhs, DenseShape, ScaledPermutationShape, ProductTag>;
 
  template <typename Dest>
  static EIGEN_DEVICE_FUNC void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) {
    const DenseType dense(lhs);
    DenseImpl::evalTo(dst, dense, rhs);
  }
  template <typename Dest>
  static EIGEN_DEVICE_FUNC void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) {
    const DenseType dense(lhs);
    DenseImpl::scaleAndAddTo(dst, dense, rhs, alpha);
  }
};
 
template <typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Rhs, ScaledPermutationShape, TriangularShape, ProductTag>
    : scaled_permutation_view_product_impl<Lhs, Rhs, ProductTag, true> {};
template <typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Rhs, ScaledPermutationShape, SelfAdjointShape, ProductTag>
    : scaled_permutation_view_product_impl<Lhs, Rhs, ProductTag, true> {};
template <typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Rhs, TriangularShape, ScaledPermutationShape, ProductTag>
    : scaled_permutation_view_product_impl<Lhs, Rhs, ProductTag, false> {};
template <typename Lhs, typename Rhs, int ProductTag>
struct generic_product_impl<Lhs, Rhs, SelfAdjointShape, ScaledPermutationShape, ProductTag>
    : scaled_permutation_view_product_impl<Lhs, Rhs, ProductTag, false> {};
 
}  // namespace internal
 
}  // namespace Eigen
 
#endif  // EIGEN_SCALEDPERMUTATIONMATRIX_H