Transpositions.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010-2011 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/.
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_TRANSPOSITIONS_H
#define EIGEN_TRANSPOSITIONS_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
template <typename Derived>
class TranspositionsBase {
  using Traits = internal::traits<Derived>;
 
 public:
  using IndicesType = typename Traits::IndicesType;
  using StorageIndex = typename IndicesType::Scalar;
  using Index = Eigen::Index;  
 
  EIGEN_DEVICE_FUNC Derived& derived() { return *static_cast<Derived*>(this); }
  EIGEN_DEVICE_FUNC const Derived& derived() const { return *static_cast<const Derived*>(this); }

  template <typename OtherDerived>
  Derived& operator=(const TranspositionsBase<OtherDerived>& other) {
    indices() = other.indices();
    return derived();
  }

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

  EIGEN_DEVICE_FUNC inline const StorageIndex& coeff(Index i) const { return indices().coeff(i); }
  inline StorageIndex& coeffRef(Index i) { return indices().coeffRef(i); }
  inline const StorageIndex& operator()(Index i) const { return indices()(i); }
  inline StorageIndex& operator()(Index i) { return indices()(i); }
  inline const StorageIndex& operator[](Index i) const { return indices()(i); }
  inline StorageIndex& operator[](Index i) { return indices()(i); }

  EIGEN_DEVICE_FUNC const IndicesType& indices() const { return derived().indices(); }
  EIGEN_DEVICE_FUNC IndicesType& indices() { return derived().indices(); }

  inline void resize(Index newSize) { indices().resize(newSize); }

  void setIdentity() {
    for (StorageIndex i = 0; i < indices().size(); ++i) coeffRef(i) = i;
  }
 
  // FIXME: do we want such methods ?
  // might be useful when the target matrix expression is complex, e.g.:
  // object.matrix().block(..,..,..,..) = trans * object.matrix().block(..,..,..,..);
  /*
  template<typename MatrixType>
  void applyForwardToRows(MatrixType& mat) const
  {
    for(Index k=0 ; k<size() ; ++k)
      if(m_indices(k)!=k)
        mat.row(k).swap(mat.row(m_indices(k)));
  }
 
  template<typename MatrixType>
  void applyBackwardToRows(MatrixType& mat) const
  {
    for(Index k=size()-1 ; k>=0 ; --k)
      if(m_indices(k)!=k)
        mat.row(k).swap(mat.row(m_indices(k)));
  }
  */

  inline Transpose<TranspositionsBase> inverse() const { return Transpose<TranspositionsBase>(derived()); }

  inline Transpose<TranspositionsBase> transpose() const { return Transpose<TranspositionsBase>(derived()); }

  Transpose<TranspositionsBase> adjoint() const { return Transpose<TranspositionsBase>(derived()); }
};
 
namespace internal {
template <int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_>
struct traits<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_> >
    : traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_> > {
  using IndicesType = Matrix<StorageIndex_, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>;
  using StorageKind = TranspositionsStorage;
};
}  // namespace internal

 
template <int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_>
class Transpositions
    : public TranspositionsBase<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_> > {
  using Traits = internal::traits<Transpositions>;
 
 public:
  using Base = TranspositionsBase<Transpositions>;
  using IndicesType = typename Traits::IndicesType;
  using StorageIndex = typename IndicesType::Scalar;
 
  inline Transpositions() {}

  template <typename OtherDerived>
  inline Transpositions(const TranspositionsBase<OtherDerived>& other) : m_indices(other.indices()) {}

  template <typename Other>
  explicit inline Transpositions(const MatrixBase<Other>& indices) : m_indices(indices) {}

  template <typename OtherDerived>
  Transpositions& operator=(const TranspositionsBase<OtherDerived>& other) {
    return Base::operator=(other);
  }

  inline Transpositions(Index size) : m_indices(size) {}

  EIGEN_DEVICE_FUNC const IndicesType& indices() const { return m_indices; }
  EIGEN_DEVICE_FUNC IndicesType& indices() { return m_indices; }
 
 protected:
  IndicesType m_indices;
};
 
namespace internal {
template <int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_, int PacketAccess_>
struct traits<Map<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>, PacketAccess_> >
    : traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_> > {
  using IndicesType = Map<const Matrix<StorageIndex_, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, PacketAccess_>;
  using StorageIndex = StorageIndex_;
  using StorageKind = TranspositionsStorage;
};
}  // namespace internal
 
template <int SizeAtCompileTime, int MaxSizeAtCompileTime, typename StorageIndex_, int PacketAccess>
class Map<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>, PacketAccess>
    : public TranspositionsBase<
          Map<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime, StorageIndex_>, PacketAccess> > {
  using Traits = internal::traits<Map>;
 
 public:
  using Base = TranspositionsBase<Map>;
  using IndicesType = typename Traits::IndicesType;
  using StorageIndex = typename IndicesType::Scalar;
 
  explicit inline Map(const StorageIndex* indicesPtr) : m_indices(indicesPtr) {}
 
  inline Map(const StorageIndex* indicesPtr, Index size) : m_indices(indicesPtr, size) {}

  template <typename OtherDerived>
  Map& operator=(const TranspositionsBase<OtherDerived>& other) {
    return Base::operator=(other);
  }
 
#ifndef EIGEN_PARSED_BY_DOXYGEN
  Map& operator=(const Map& other) {
    m_indices = other.m_indices;
    return *this;
  }
#endif

  EIGEN_DEVICE_FUNC const IndicesType& indices() const { return m_indices; }

  EIGEN_DEVICE_FUNC IndicesType& indices() { return m_indices; }
 
 protected:
  IndicesType m_indices;
};
 
namespace internal {
template <typename IndicesType_>
struct traits<TranspositionsWrapper<IndicesType_> > : traits<PermutationWrapper<IndicesType_> > {
  using StorageKind = TranspositionsStorage;
};
}  // namespace internal
 
template <typename IndicesType_>
class TranspositionsWrapper : public TranspositionsBase<TranspositionsWrapper<IndicesType_> > {
  using Traits = internal::traits<TranspositionsWrapper>;
 
 public:
  using Base = TranspositionsBase<TranspositionsWrapper>;
  using IndicesType = typename Traits::IndicesType;
  using StorageIndex = typename IndicesType::Scalar;
 
  explicit inline TranspositionsWrapper(IndicesType& indices) : m_indices(indices) {}

  template <typename OtherDerived>
  TranspositionsWrapper& operator=(const TranspositionsBase<OtherDerived>& other) {
    return Base::operator=(other);
  }

  EIGEN_DEVICE_FUNC const IndicesType& indices() const { return m_indices; }

  EIGEN_DEVICE_FUNC IndicesType& indices() { return m_indices; }
 
 protected:
  typename IndicesType::Nested m_indices;
};

template <typename MatrixDerived, typename TranspositionsDerived>
EIGEN_DEVICE_FUNC const Product<MatrixDerived, TranspositionsDerived, AliasFreeProduct> operator*(
    const MatrixBase<MatrixDerived>& matrix, const TranspositionsBase<TranspositionsDerived>& transpositions) {
  return Product<MatrixDerived, TranspositionsDerived, AliasFreeProduct>(matrix.derived(), transpositions.derived());
}

template <typename TranspositionsDerived, typename MatrixDerived>
EIGEN_DEVICE_FUNC const Product<TranspositionsDerived, MatrixDerived, AliasFreeProduct> operator*(
    const TranspositionsBase<TranspositionsDerived>& transpositions, const MatrixBase<MatrixDerived>& matrix) {
  return Product<TranspositionsDerived, MatrixDerived, AliasFreeProduct>(transpositions.derived(), matrix.derived());
}
 
// Template partial specialization for transposed/inverse transpositions
 
namespace internal {
 
template <typename Derived>
struct traits<Transpose<TranspositionsBase<Derived> > > : traits<Derived> {};
 
}  // end namespace internal
 
template <typename TranspositionsDerived>
class Transpose<TranspositionsBase<TranspositionsDerived> > {
  using TranspositionType = TranspositionsDerived;
  using IndicesType = typename TranspositionType::IndicesType;
 
 public:
  explicit Transpose(const TranspositionType& t) : m_transpositions(t) {}
 
  EIGEN_DEVICE_FUNC constexpr Index size() const noexcept { return m_transpositions.size(); }
  EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return m_transpositions.size(); }
  EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return m_transpositions.size(); }

  template <typename OtherDerived>
  friend const Product<OtherDerived, Transpose, AliasFreeProduct> operator*(const MatrixBase<OtherDerived>& matrix,
                                                                            const Transpose& trt) {
    return Product<OtherDerived, Transpose, AliasFreeProduct>(matrix.derived(), trt);
  }

  template <typename OtherDerived>
  const Product<Transpose, OtherDerived, AliasFreeProduct> operator*(const MatrixBase<OtherDerived>& matrix) const {
    return Product<Transpose, OtherDerived, AliasFreeProduct>(*this, matrix.derived());
  }
 
  EIGEN_DEVICE_FUNC const TranspositionType& nestedExpression() const { return m_transpositions; }
 
 protected:
  const TranspositionType& m_transpositions;
};
 
}  // end namespace Eigen
 
#endif  // EIGEN_TRANSPOSITIONS_H