StructuredMatrixUtils.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
 
// References:
//  [1] N. J. Higham, "Accuracy and Stability of Numerical Algorithms", 2nd ed.,
//      SIAM, 2002, chapter 27. Avoiding spurious overflow by rescaling with
//      powers of two, the technique behind structured_exponent_bound() and the
//      balanced determinant accumulations.
//  [2] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974.
//      Scaling by a power of two is exact, so the balanced accumulations
//      introduce no roundoff of their own.
 
#ifndef EIGEN_STRUCTURED_MATRIX_UTILS_H
#define EIGEN_STRUCTURED_MATRIX_UTILS_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
namespace internal {
 
// A distinct shape routes assignment through evalTo/addTo/subTo and lets one
// product specialization cover every product tag without colliding with DenseShape.
struct StructuredShape {};
 
// Below this dimension the FFT setup costs more than a plain O(n^2) evaluation,
// so the structured operators fall back to a direct segment-based product.
constexpr Index structured_direct_threshold() { return 32; }
 
// Below this dimension even the segment-based direct product loses to a plain
// scalar loop: the per-segment setup dominates when the average segment holds
// fewer than a couple of packets (measured crossover on AVX2 hardware).
constexpr Index structured_scalar_threshold() { return 16; }
 
// Numerical scale exponents are independent of Eigen's configurable dimension
// index. A 32-bit Index can overflow while accumulating O(n^2) factor
// exponents for dimensions that are otherwise practical.
using structured_exponent_type = numext::int64_t;

template <typename Scalar>
typename NumTraits<Scalar>::Real structured_component_magnitude_impl(const Scalar& z, std::true_type) {
  using RealScalar = typename NumTraits<Scalar>::Real;
  using Binary = binary_floating_point_traits<RealScalar>;
  return numext::bit_cast<RealScalar>(
      larger_magnitude_bits<RealScalar>(Binary::magnitude(numext::real(z)), Binary::magnitude(numext::imag(z))));
}
template <typename Scalar>
typename NumTraits<Scalar>::Real structured_component_magnitude_impl(const Scalar& z, std::false_type) {
  return numext::maxi(numext::abs(numext::real(z)), numext::abs(numext::imag(z)));
}
template <typename Scalar>
typename NumTraits<Scalar>::Real structured_component_magnitude(const Scalar& z) {
  using RealScalar = typename NumTraits<Scalar>::Real;
  return structured_component_magnitude_impl(z,
                                             bool_constant < complex_array_access<Scalar>::value &&
                                                 use_subnormal_preserving_scaling<RealScalar, RealScalar>::value > ());
}

template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct structured_balance_impl {
  using RealScalar = typename NumTraits<Scalar>::Real;
  template <typename Exponent>
  static Scalar run(const Scalar& z, Exponent& exponent) {
    const RealScalar mag = structured_component_magnitude(z);
    if (numext::is_exactly_zero_no_flush(mag) || !(numext::isfinite)(mag)) return z;
    const int e = frexp_exponent_preserving_subnormals(mag);
    exponent += e;
    return apply_exponent(z, -e);
  }
  static Scalar apply_exponent(const Scalar& z, int e) {
    return Scalar(ldexp_preserving_subnormals(numext::real(z), e), ldexp_preserving_subnormals(numext::imag(z), e));
  }
};
 
template <typename Scalar>
struct structured_balance_impl<Scalar, false> {
  template <typename Exponent>
  static Scalar run(const Scalar& x, Exponent& exponent) {
    if (numext::is_exactly_zero_no_flush(x) || !(numext::isfinite)(x)) return x;
    const int e = frexp_exponent_preserving_subnormals(x);
    exponent += e;
    return apply_exponent(x, -e);
  }
  static Scalar apply_exponent(const Scalar& x, int e) { return ldexp_preserving_subnormals(x, e); }
};
 
template <typename Scalar, typename Exponent>
Scalar structured_balance(const Scalar& z, Exponent& exponent) {
  return structured_balance_impl<Scalar>::run(z, exponent);
}

template <typename Scalar, typename Exponent>
Scalar structured_ldexp_clamped(const Scalar& z, Exponent exponent) {
  constexpr Exponent kMaxExponent = Exponent(1) << 24;
  const int e = static_cast<int>(numext::mini(numext::maxi(exponent, -kMaxExponent), kMaxExponent));
  return structured_balance_impl<Scalar>::apply_exponent(z, e);
}

template <typename RealVectorType>
std::vector<Index> structured_svd_permutation(const RealVectorType& mods) {
  using RealScalar = typename RealVectorType::Scalar;
  std::vector<Index> perm;
  perm.reserve(static_cast<std::size_t>(mods.size()));
  for (Index k = 0; k < mods.size(); ++k) perm.push_back(k);
  std::stable_sort(perm.begin(), perm.end(), [&mods](Index a, Index b) {
    const RealScalar ka = mods[a], kb = mods[b];
    // isgreater is quiet for NaNs; the second clause implements NaN-last ordering.
    return std::isgreater(ka, kb) || (!(numext::isnan)(ka) && (numext::isnan)(kb));
  });
  return perm;
}

#ifndef EIGEN_AVOID_THREAD_LOCAL
template <typename RealScalar>
FFT<RealScalar>& structured_fft_engine() {
  static thread_local FFT<RealScalar> fft;
  return fft;
}
#else
template <typename RealScalar>
FFT<RealScalar> structured_fft_engine() {
  return FFT<RealScalar>();
}
#endif

inline Index fft_next_good_size(Index n) {
  if (n < 1) return 1;
  for (Index m = n;; ++m) {
    Index r = m;
    while (r % 2 == 0) r /= 2;
    while (r % 3 == 0) r /= 3;
    while (r % 5 == 0) r /= 5;
    if (r == 1) return m;
  }
}

template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
struct structured_scalar_part_impl {
  template <typename Xpr>
  static const Xpr& run(const Xpr& xpr) {
    return xpr;
  }
  static const Scalar& run_scalar(const Scalar& x) { return x; }
};
 
template <typename Scalar>
struct structured_scalar_part_impl<Scalar, false> {
  template <typename Xpr>
  static typename Xpr::RealReturnType run(const Xpr& xpr) {
    return xpr.real();
  }
  static Scalar run_scalar(const std::complex<Scalar>& x) { return numext::real(x); }
};

template <typename Xpr>
bool structured_exponent_bound_finite(const Xpr& x, int& e) {
  using ScalarTraits = NumTraits<typename Xpr::Scalar>;
  using RealScalar = typename ScalarTraits::Real;
  // maxCoeff() asserts on an empty input, and a degenerate operand -- a rank-0
  // factor, a solve with no right-hand sides -- reaches here legitimately.
  // An empty operand bounds nothing, so its exponent bound is 0.
  e = 0;
  if (x.size() == 0) return true;
  RealScalar m;
  if (ScalarTraits::IsComplex)
    // realView() reduces over both components in one pass, vectorized for
    // direct-access storage; the strided real()/imag() views never vectorize.
    m = x.realView().cwiseAbs().maxCoeff();
  else
    m = x.cwiseAbs().maxCoeff();
  // A SIMD unit that flushes subnormal inputs (ARMv7 NEON, Arm FZ, DAZ) reads
  // an all-subnormal operand as zero; the rescan recovers the largest component
  // from its representation.
  m = safe_scaling<RealScalar>::recover_flushed_max_coeff(x, m);
  e = 0;
  if (!(numext::isfinite)(m)) return false;
  if (!numext::is_exactly_zero_no_flush(m)) {
    e = frexp_exponent_preserving_subnormals(m);
    if (ScalarTraits::IsComplex) ++e;
  }
  return true;
}

template <typename Xpr>
int structured_exponent_bound(const Xpr& x) {
  int e;
  structured_exponent_bound_finite(x, e);
  return e;
}
 
// The packet form of M *= 2^e. ldexp saturates entrywise for every e without
// forming a possibly unrepresentable 2^e; for std::complex storage, realView()
// exposes both components to the real ldexp packets.
template <typename Xpr, std::enable_if_t<!NumTraits<typename Xpr::Scalar>::IsComplex, bool> = true>
void structured_ldexp_entries_packet(Xpr& M, int e) {
  M.array() = M.array().ldexp(e);
}
 
template <typename Xpr, std::enable_if_t<complex_array_access<typename Xpr::Scalar>::value, bool> = true>
void structured_ldexp_entries_packet(Xpr& M, int e) {
  M.realView().array() = M.realView().array().ldexp(e);
}
 
template <typename Xpr, std::enable_if_t<NumTraits<typename Xpr::Scalar>::IsComplex &&
                                             !complex_array_access<typename Xpr::Scalar>::value,
                                         bool> = true>
void structured_ldexp_entries_packet(Xpr& M, int e) {
  using Scalar = typename Xpr::Scalar;
  M = M.unaryExpr([e](const Scalar& z) { return structured_ldexp_clamped(z, Index(e)); });
}
 
template <typename Xpr>
void structured_ldexp_entries_impl(Xpr& M, int e, int, std::false_type) {
  structured_ldexp_entries_packet(M, e);
}
 
template <typename Xpr>
void structured_ldexp_entries_impl(Xpr& M, int e, int bound, std::true_type) {
  using Scalar = typename Xpr::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  // The largest component is at least 2^(componentBound - 1): a complex bound
  // carries one extra bit for the modulus. The input holds no significant
  // subnormal when componentBound - 1 >= recovery, and neither does the result
  // when componentBound - 1 + e >= recovery.
  const int componentBound = bound - (NumTraits<Scalar>::IsComplex ? 1 : 0);
  const int recovery = safe_scaling<RealScalar>::subnormal_recovery_exponent();
  if (componentBound - 1 >= recovery && componentBound - 1 + e >= recovery)
    structured_ldexp_entries_packet(M, e);
  else
    M = M.unaryExpr(scale_by_exponent_op<RealScalar>(e));
}

template <typename Xpr>
void structured_ldexp_entries(Xpr& M, int e, int bound) {
  if (e == 0) return;
  using Scalar = typename Xpr::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  structured_ldexp_entries_impl(M, e, bound,
                                bool_constant<use_subnormal_preserving_scaling<RealScalar, Scalar>::value>());
}

template <typename Xpr>
void structured_ldexp_entries(Xpr& M, int e) {
  if (e == 0) return;
  structured_ldexp_entries(M, e, structured_exponent_bound(M));
}
 
template <typename Xpr>
void structured_ldexp_entries_exact_impl(Xpr& M, int e, std::false_type) {
  structured_ldexp_entries_packet(M, e);
}
template <typename Xpr>
void structured_ldexp_entries_exact_impl(Xpr& M, int e, std::true_type) {
  using RealScalar = typename NumTraits<typename Xpr::Scalar>::Real;
  M = M.unaryExpr(scale_by_exponent_op<RealScalar>(e));
}

template <typename Xpr>
void structured_ldexp_entries_exact(Xpr& M, int e) {
  if (e == 0) return;
  using Scalar = typename Xpr::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  structured_ldexp_entries_exact_impl(M, e,
                                      bool_constant<use_subnormal_preserving_scaling<RealScalar, Scalar>::value>());
}

template <typename ComplexVectorType>
ComplexVectorType structured_reverse_symbol(const ComplexVectorType& symbol) {
  const Index p = symbol.size();
  ComplexVectorType reversed(p);
  if (p > 0) {
    reversed[0] = symbol[0];
    reversed.tail(p - 1) = symbol.tail(p - 1).reverse();
  }
  return reversed;
}

template <typename Scalar, typename Dest, typename Rhs>
void structured_fft_apply(Dest& dst, const Matrix<std::complex<typename NumTraits<Scalar>::Real>, Dynamic, 1>& symbol,
                          Index outSize, const Rhs& rhs, const Scalar& alpha) {
  using RealScalar = typename NumTraits<Scalar>::Real;
  using Complex = std::complex<RealScalar>;
  using ComplexVector = Matrix<Complex, Dynamic, 1>;
 
  const Index p = symbol.size();
  eigen_assert(rhs.rows() <= p && outSize <= p);
  if (p == 1) {
    // The length-one DFT is the identity and is unsupported by kissfft.
    dst.row(0) += alpha * structured_scalar_part_impl<Scalar>::run(Complex(symbol.coeff(0)) *
                                                                   rhs.row(0).template cast<Complex>());
    return;
  }
 
  auto&& fft = structured_fft_engine<RealScalar>();
  ComplexVector xt = ComplexVector::Zero(p);
  ComplexVector xf(p), yt(p);
  for (Index k = 0; k < rhs.cols(); ++k) {
    xt.head(rhs.rows()) = rhs.col(k).template cast<Complex>();
    fft.fwd(xf, xt, p);
    xf.array() *= symbol.array();
    fft.inv(yt, xf, p);
    dst.col(k) += alpha * structured_scalar_part_impl<Scalar>::run(yt.head(outSize));
  }
}

template <typename SymbolType, typename ModsType, typename RealScalar>
SymbolType structured_pinv_symbol(const SymbolType& symbol, const ModsType& mods, const RealScalar& tol) {
  using Complex = typename SymbolType::Scalar;
  const auto w = (mods.array() < tol).select(RealScalar(0), mods.array().inverse()).template cast<Complex>().eval();
  return (symbol.array().conjugate() * w * w).matrix();
}

template <typename SymbolType, typename ModsType>
typename ModsType::Scalar structured_rank_threshold(const SymbolType& symbol, const ModsType& mods) {
  using RealScalar = typename ModsType::Scalar;
  const RealScalar factor = RealScalar(mods.size()) * NumTraits<RealScalar>::epsilon();
  RealScalar tol = factor * mods.maxCoeff();
  if (!(numext::isfinite)(tol)) tol = (RealScalar(2) * factor) * (symbol * RealScalar(0.5)).cwiseAbs().maxCoeff();
  return numext::maxi(tol, (std::numeric_limits<RealScalar>::min)());
}

template <typename Op, typename Rhs>
struct structured_product_impl : generic_product_impl_base<Op, Rhs, structured_product_impl<Op, Rhs>> {
  using Scalar = typename Product<Op, Rhs>::Scalar;
 
  template <typename Dest>
  static void evalTo(Dest& dst, const Op& lhs, const Rhs& rhs) {
    dst.setZero();
    scaleAndAddTo(dst, lhs, rhs, Scalar(1));
  }
 
  template <typename Dest>
  static void scaleAndAddTo(Dest& dst, const Op& lhs, const Rhs& rhs, const Scalar& alpha) {
    using RhsNested = typename nested_eval<Rhs, Op::RowsAtCompileTime>::type;
    RhsNested actualRhs(rhs);
    lhs.addProduct(dst, actualRhs, alpha);
  }
};
 
}  // namespace internal
 
}  // namespace Eigen
 
#endif  // EIGEN_STRUCTURED_MATRIX_UTILS_H