KroneckerOperator.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] C. F. Van Loan, "The ubiquitous Kronecker product", Journal of
//      Computational and Applied Mathematics 123 (2000), 85-100. The vec
//      identity, with vec stacking columns, driving every product and solve
//      below -- stated there as
//      Y = C X B^T <=> vec(Y) = (B (x) C) vec(X), i.e. with this file's
//      operand naming (A (x) B) vec(X) = vec(B X A^T) -- together with the
//      factor-wise identities for the inverse, pseudo-inverse,
//      eigendecomposition, SVD and determinant of a Kronecker product.
//  [2] 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 the exponent-balanced determinant.
//  [3] P. H. Sterbenz, "Floating-Point Computation", Prentice-Hall, 1974.
//      Scaling by a power of two is exact, the property the determinant's
//      rescaling relies on.
 
#ifndef EIGEN_STRUCTURED_KRONECKER_OPERATOR_H
#define EIGEN_STRUCTURED_KRONECKER_OPERATOR_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
template <typename LhsMatrix, typename RhsMatrix>
class KroneckerOperator;
template <typename LhsMatrix, typename RhsMatrix>
class KroneckerSum;
 
namespace internal {
 
template <typename LhsMatrix, typename RhsMatrix>
struct traits<KroneckerOperator<LhsMatrix, RhsMatrix>> {
  using Scalar = typename LhsMatrix::Scalar;
  using StorageKind = Dense;
  using XprKind = MatrixXpr;
  using StorageIndex = int;
  // size_at_compile_time is the compile-time dimension product, Dynamic when a
  // factor is unknown or the product would overflow int.
  static constexpr int RowsAtCompileTime =
      size_at_compile_time(traits<LhsMatrix>::RowsAtCompileTime, traits<RhsMatrix>::RowsAtCompileTime);
  static constexpr int ColsAtCompileTime =
      size_at_compile_time(traits<LhsMatrix>::ColsAtCompileTime, traits<RhsMatrix>::ColsAtCompileTime);
  static constexpr int MaxRowsAtCompileTime = RowsAtCompileTime;
  static constexpr int MaxColsAtCompileTime = ColsAtCompileTime;
  // Deliberately no NestByRefBit: transpose(), conjugate(), adjoint(), inverse(),
  // eigenvectors(), matrixU() and matrixV() return owning temporaries (the
  // operator stores its factors by value), so Product must nest the operator by
  // value for a delayed-evaluated product expression to keep its left factor
  // alive. The copy is O(m1 n1 + m2 n2), negligible against the product
  // evaluation.
  static constexpr int Flags = 0;
};
 
template <typename LhsMatrix, typename RhsMatrix>
struct evaluator_traits<KroneckerOperator<LhsMatrix, RhsMatrix>> {
  using Kind = IndexBased;
  using Shape = StructuredShape;
};
 
// Factor kinds, each with its operations in the kron_factor_* helpers below:
// dense; diagonal (stored as its diagonal); sparse (compressed, solved by
// SparseLU); identity (a kron_identity_factor, dimensions only, skipped in
// products); Kronecker (a nested KroneckerOperator, for three or more factors);
// Kronecker sum (a KroneckerSum, see KroneckerSum.h).
 
template <typename Scalar_, int Rows_, int Cols_>
class kron_identity_factor;
 
template <typename Scalar_, int Rows_, int Cols_>
struct traits<kron_identity_factor<Scalar_, Rows_, Cols_>> : traits<Matrix<Scalar_, Rows_, Cols_>> {};

template <typename Scalar_, int Rows_, int Cols_>
class kron_identity_factor : public EigenBase<kron_identity_factor<Scalar_, Rows_, Cols_>> {
 public:
  using Scalar = Scalar_;
  using PlainObject = Matrix<Scalar, Rows_, Cols_>;
  static constexpr int RowsAtCompileTime = Rows_;
  static constexpr int ColsAtCompileTime = Cols_;
 
  kron_identity_factor(Index rows, Index cols) : m_rows(rows), m_cols(cols) {}
  template <typename PlainObjectType>
  explicit kron_identity_factor(const CwiseNullaryOp<scalar_identity_op<Scalar>, PlainObjectType>& identity)
      : m_rows(identity.rows()), m_cols(identity.cols()) {}
 
  EIGEN_DEVICE_FUNC constexpr Index rows() const { return m_rows.value(); }
  EIGEN_DEVICE_FUNC constexpr Index cols() const { return m_cols.value(); }
 
 private:
  variable_if_dynamic<Index, Rows_> m_rows;
  variable_if_dynamic<Index, Cols_> m_cols;
};
 
template <typename Factor>
struct kron_factor_is_diagonal : std::false_type {};
template <typename Scalar, int Size, int MaxSize>
struct kron_factor_is_diagonal<DiagonalMatrix<Scalar, Size, MaxSize>> : std::true_type {};
 
template <typename Factor>
struct kron_factor_is_dense_matrix : std::false_type {};
template <typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
struct kron_factor_is_dense_matrix<Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols>> : std::true_type {};
 
template <typename Factor>
struct kron_factor_is_sparse_matrix : std::false_type {};
template <typename Scalar, int Options, typename StorageIndex>
struct kron_factor_is_sparse_matrix<SparseMatrix<Scalar, Options, StorageIndex>> : std::true_type {};
 
template <typename Factor>
struct kron_factor_is_identity : std::false_type {};
template <typename Scalar, int Rows, int Cols>
struct kron_factor_is_identity<kron_identity_factor<Scalar, Rows, Cols>> : std::true_type {};
 
template <typename Factor>
struct kron_factor_is_kronecker : std::false_type {};
template <typename LhsMatrix, typename RhsMatrix>
struct kron_factor_is_kronecker<KroneckerOperator<LhsMatrix, RhsMatrix>> : std::true_type {};
 
template <typename Factor>
struct kron_factor_is_kronecker_sum : std::false_type {};
template <typename LhsMatrix, typename RhsMatrix>
struct kron_factor_is_kronecker_sum<KroneckerSum<LhsMatrix, RhsMatrix>> : std::true_type {};
 
// The factor kind, the dispatch key of kron_factor_ops and kron_factor_solver.
constexpr int kKronDenseFactor = 0;
constexpr int kKronDiagonalFactor = 1;
constexpr int kKronSparseFactor = 2;
constexpr int kKronIdentityFactor = 3;
constexpr int kKronKroneckerFactor = 4;
constexpr int kKronSumFactor = 5;
 
template <typename Factor>
constexpr int kron_factor_kind() {
  return kron_factor_is_diagonal<Factor>::value        ? kKronDiagonalFactor
         : kron_factor_is_sparse_matrix<Factor>::value ? kKronSparseFactor
         : kron_factor_is_identity<Factor>::value      ? kKronIdentityFactor
         : kron_factor_is_kronecker<Factor>::value     ? kKronKroneckerFactor
         : kron_factor_is_kronecker_sum<Factor>::value ? kKronSumFactor
                                                       : kKronDenseFactor;
}

template <typename Derived, bool StoredAsIs = kron_factor_is_identity<Derived>::value ||
                                              kron_factor_is_kronecker<Derived>::value ||
                                              kron_factor_is_kronecker_sum<Derived>::value>
struct kron_factor_storage {
  using type = typename Derived::PlainObject;
};
template <typename Derived>
struct kron_factor_storage<Derived, true> {
  using type = Derived;
};
template <typename Scalar, typename PlainObjectType>
struct kron_factor_storage<CwiseNullaryOp<scalar_identity_op<Scalar>, PlainObjectType>, false> {
  using type = kron_identity_factor<Scalar, PlainObjectType::RowsAtCompileTime, PlainObjectType::ColsAtCompileTime>;
};

template <typename Factor, int Kind = kron_factor_kind<Factor>()>
struct kron_factor_visitable {
  using type = Factor;
  static const Factor& get(const Factor& f) { return f; }
};

template <typename LhsMatrix, typename RhsMatrix,
          bool Unchanged = std::is_same<typename kron_factor_visitable<LhsMatrix>::type, LhsMatrix>::value &&
                           std::is_same<typename kron_factor_visitable<RhsMatrix>::type, RhsMatrix>::value>
struct kron_operator_visitable {
  using Factor = KroneckerOperator<LhsMatrix, RhsMatrix>;
  using type = Factor;
  static const Factor& get(const Factor& f) { return f; }
};
template <typename LhsMatrix, typename RhsMatrix>
struct kron_operator_visitable<LhsMatrix, RhsMatrix, false> {
  using LhsVisitable = kron_factor_visitable<LhsMatrix>;
  using RhsVisitable = kron_factor_visitable<RhsMatrix>;
  using type = KroneckerOperator<typename LhsVisitable::type, typename RhsVisitable::type>;
  static type get(const KroneckerOperator<LhsMatrix, RhsMatrix>& f) {
    return type(LhsVisitable::get(f.lhs()), RhsVisitable::get(f.rhs()));
  }
};
template <typename LhsMatrix, typename RhsMatrix>
struct kron_factor_visitable<KroneckerOperator<LhsMatrix, RhsMatrix>, kKronKroneckerFactor>
    : kron_operator_visitable<LhsMatrix, RhsMatrix> {};

template <typename Stacked, typename Xpr>
void kron_stack_columns(Stacked& stacked, const Xpr& x, Index rows, Index cols) {
  stacked.resize(rows * x.cols(), cols);
  for (Index k = 0; k < x.cols(); ++k)
    stacked.middleRows(k * rows, rows) = x.col(k).reshaped(rows, cols).template cast<typename Stacked::Scalar>();
}

template <typename WorkScalar>
Index kron_rhs_chunk(Index perColumn) {
  const Index budget = Index(l2CacheSize() / 8) / Index(sizeof(WorkScalar));
  return numext::maxi(Index(1), budget / numext::maxi(perColumn, Index(1)));
}
 
template <typename Factor, int Kind = kron_factor_kind<Factor>()>
struct kron_factor_ops {
  // Dense factor.
  using Scalar = typename Factor::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using TransposedFactor = Matrix<Scalar, Factor::ColsAtCompileTime, Factor::RowsAtCompileTime>;
  using InverseFactor = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
  // Every entry is stored: a materialization has no structurally zero blocks
  // to clear, and forEachNonZero visits all m*n entries in column-major order.
  static constexpr bool StoresAllEntries = true;
 
  static void prepare(Factor&) {}
  static bool coeffIfStored(const Factor& f, Index row, Index col, Scalar& value) {
    value = f.coeff(row, col);
    return true;
  }
  template <typename Visitor>
  static void forEachNonZero(const Factor& f, Visitor&& visit) {
    for (Index j = 0; j < f.cols(); ++j)
      for (Index i = 0; i < f.rows(); ++i) visit(i, j, f.coeff(i, j));
  }
  static Matrix<Index, Dynamic, 1> innerNonZeros(const Factor& f, bool rowMajor) {
    return Matrix<Index, Dynamic, 1>::Constant(rowMajor ? f.rows() : f.cols(), rowMajor ? f.cols() : f.rows());
  }
  static auto transposed(const Factor& f) { return f.transpose(); }
  static auto conjugated(const Factor& f) { return f.conjugate(); }
  static auto adjointed(const Factor& f) { return f.adjoint(); }
  static auto inversed(const Factor& f) { return f.inverse(); }
  // The factor as a dense expression, for the decomposition family.
  static const Factor& denseFactor(const Factor& f) { return f; }
  // The factor as the operand of a dense block assignment dst = a * f.
  static const Factor& blockOperand(const Factor& f) { return f; }
  static bool isSquareIdentity(const Factor&) { return false; }
  // dst += alpha F X and dst += alpha X F^T: the two sides of the vec-trick product B X A^T.
  // work: scratch for a nested factor's right product, owned by the caller's loop.
  template <typename Dst, typename Alpha, typename Xpr>
  static void addLeftProduct(Dst& dst, const Alpha& alpha, const Factor& f, const Xpr& X) {
    dst.noalias() += alpha * (f * X);
  }
  template <typename Dst, typename Alpha, typename Xpr, typename Work>
  static void addRightProduct(Dst& dst, const Alpha& alpha, const Xpr& X, const Factor& f, Work&) {
    dst.noalias() += alpha * (X * f.transpose());
  }
  static int exponentBound(const Factor& f) { return structured_exponent_bound(f); }
  static Scalar balancedDet(const Factor& M, Index& exponent) {
    // Scale small factors up, exactly, so that the elimination runs on normal
    // numbers: a subnormal near 2^e carries only e - min_exponent + digits bits.
    // Scaling down could erase small pivots in factors with a wide exponent range,
    // and every entry is scaled exactly because a small entry can be a pivot.
    const int bound = exponentBound(M);
    const int scaleExponent = numext::mini(bound, 0);
    Factor normalized = M;
    structured_ldexp_entries_exact(normalized, -scaleExponent);
    exponent += M.rows() * scaleExponent;
    PartialPivLU<Factor> lu(normalized);
    Scalar m = Scalar(RealScalar(lu.permutationP().determinant()));  // +-1
    for (Index i = 0; i < M.rows(); ++i)
      m = structured_balance(m * structured_balance(lu.matrixLU().coeff(i, i), exponent), exponent);
    return m;
  }
};
 
template <typename Factor>
struct kron_factor_ops<Factor, kKronDiagonalFactor> {
  // Diagonal factor: everything runs on the stored diagonal vector.
  using Scalar = typename Factor::Scalar;
  using TransposedFactor = typename Factor::PlainObject;  // a diagonal matrix is its own transpose
  using InverseFactor = typename Factor::PlainObject;     // entrywise reciprocals stay diagonal
  static constexpr bool StoresAllEntries = false;
 
  static void prepare(Factor&) {}
  static bool coeffIfStored(const Factor& f, Index row, Index col, Scalar& value) {
    if (row != col) return false;
    value = f.diagonal().coeff(row);
    return true;
  }
  template <typename Visitor>
  static void forEachNonZero(const Factor& f, Visitor&& visit) {
    for (Index k = 0; k < f.rows(); ++k) visit(k, k, f.diagonal().coeff(k));
  }
  static Matrix<Index, Dynamic, 1> innerNonZeros(const Factor& f, bool) {
    return Matrix<Index, Dynamic, 1>::Ones(f.rows());
  }
  static const Factor& transposed(const Factor& f) { return f; }
  static auto conjugated(const Factor& f) { return f.diagonal().conjugate().asDiagonal(); }
  static auto adjointed(const Factor& f) { return f.diagonal().conjugate().asDiagonal(); }
  static auto inversed(const Factor& f) { return f.inverse(); }
  static typename Factor::DenseMatrixType denseFactor(const Factor& f) { return f.toDenseMatrix(); }
  static const Factor& blockOperand(const Factor& f) { return f; }
  static bool isSquareIdentity(const Factor&) { return false; }
  template <typename Dst, typename Alpha, typename Xpr>
  static void addLeftProduct(Dst& dst, const Alpha& alpha, const Factor& f, const Xpr& X) {
    dst.noalias() += alpha * (f * X);
  }
  template <typename Dst, typename Alpha, typename Xpr, typename Work>
  static void addRightProduct(Dst& dst, const Alpha& alpha, const Xpr& X, const Factor& f, Work&) {
    dst.noalias() += alpha * (X * f);
  }
  static int exponentBound(const Factor& f) { return structured_exponent_bound(f.diagonal()); }
  static Scalar balancedDet(const Factor& D, Index& exponent) {
    // Scale a small diagonal up first: the balancing below compares and scales
    // in floating point, which reads a subnormal entry as zero under flush-to-zero.
    const int bound = exponentBound(D);
    const int scaleExponent = numext::mini(bound, 0);
    typename Factor::DiagonalVectorType d = D.diagonal();
    structured_ldexp_entries_exact(d, -scaleExponent);
    exponent += D.rows() * scaleExponent;
    Scalar m(1);
    for (Index i = 0; i < D.rows(); ++i) m = structured_balance(m * structured_balance(d.coeff(i), exponent), exponent);
    return m;
  }
};

template <typename SparseType>
class kron_sparse_lu : public SparseLU<SparseType> {
 public:
  using Base = SparseLU<SparseType>;
  using Scalar = typename SparseType::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
 
  Scalar balancedDet(Index& exponent) const {
    eigen_assert(Base::info() == Success);
    Scalar m = Scalar(RealScalar(Base::m_detPermR * Base::m_detPermC));
    for (Index j = 0; j < Base::cols(); ++j) {
      typename Base::SCMatrix::InnerIterator it(Base::m_Lstore, j);
      while (it && it.index() != j) ++it;
      if (it) m = structured_balance(m * structured_balance(it.value(), exponent), exponent);
    }
    return m;
  }
};
 
template <typename Factor>
struct kron_factor_ops<Factor, kKronSparseFactor> {
  // Sparse factor: the products run through the sparse-dense kernels, the
  // solve and the determinant through a SparseLU factorization.
  using Scalar = typename Factor::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using ColMajorFactor = SparseMatrix<Scalar, ColMajor, typename Factor::StorageIndex>;
  using TransposedFactor = Factor;  // the transpose is stored sparse, in the same storage order
  using InverseFactor = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;  // dense in general
  static constexpr bool StoresAllEntries = false;
 
  // coeffs() below reads the compressed value array.
  static void prepare(Factor& f) { f.makeCompressed(); }
  static bool coeffIfStored(const Factor& f, Index row, Index col, Scalar& value) {
    const Index outer = Factor::IsRowMajor ? row : col;
    const Index inner = Factor::IsRowMajor ? col : row;
    const Index start = f.outerIndexPtr()[outer], end = f.outerIndexPtr()[outer + 1];
    if (start == end) return false;
    const Index position = f.data().searchLowerIndex(start, end, inner);
    if (position == end || f.innerIndexPtr()[position] != inner) return false;
    value = f.valuePtr()[position];
    return true;
  }
  template <typename Visitor>
  static void forEachNonZero(const Factor& f, Visitor&& visit) {
    for (Index k = 0; k < f.outerSize(); ++k)
      for (typename Factor::InnerIterator it(f, k); it; ++it) visit(it.row(), it.col(), it.value());
  }
  static Matrix<Index, Dynamic, 1> innerNonZeros(const Factor& f, bool rowMajor) {
    Matrix<Index, Dynamic, 1> counts = Matrix<Index, Dynamic, 1>::Zero(rowMajor ? f.rows() : f.cols());
    for (Index k = 0; k < f.outerSize(); ++k)
      for (typename Factor::InnerIterator it(f, k); it; ++it) ++counts[rowMajor ? it.row() : it.col()];
    return counts;
  }
  static auto transposed(const Factor& f) { return f.transpose(); }
  static auto conjugated(const Factor& f) { return f.conjugate(); }
  static auto adjointed(const Factor& f) { return f.adjoint(); }
  static InverseFactor inversed(const Factor& f) {
    const ColMajorFactor colMajor(f);
    SparseLU<ColMajorFactor> lu(colMajor);
    if (lu.info() != Success) return nanMatrix(f.rows(), f.cols());
    return lu.solve(InverseFactor::Identity(f.rows(), f.cols()));
  }
  static InverseFactor nanMatrix(Index rows, Index cols) {
    return InverseFactor::Constant(rows, cols, Scalar(NumTraits<RealScalar>::quiet_NaN()));
  }
  static InverseFactor denseFactor(const Factor& f) { return InverseFactor(f); }
  static const Factor& blockOperand(const Factor& f) { return f; }
  static bool isSquareIdentity(const Factor&) { return false; }
  template <typename Dst, typename Alpha, typename Xpr>
  static void addLeftProduct(Dst& dst, const Alpha& alpha, const Factor& f, const Xpr& X) {
    dst.noalias() += alpha * (f * X);
  }
  template <typename Dst, typename Alpha, typename Xpr, typename Work>
  static void addRightProduct(Dst& dst, const Alpha& alpha, const Xpr& X, const Factor& f, Work&) {
    dst.noalias() += alpha * (X * f.transpose());
  }
  static int exponentBound(const Factor& f) { return structured_exponent_bound(f.coeffs()); }
  static Scalar balancedDet(const Factor& M, Index& exponent) {
    // Scale small factors up, exactly, so that the elimination runs on normal
    // numbers: a subnormal near 2^e carries only e - min_exponent + digits bits.
    // Scaling down could erase small pivots in factors with a wide exponent range.
    const int bound = exponentBound(M);
    const int scaleExponent = numext::mini(bound, 0);
    ColMajorFactor normalized(M);
    if (scaleExponent != 0) {
      auto values = normalized.coeffs();
      structured_ldexp_entries_exact(values, -scaleExponent);
    }
    kron_sparse_lu<ColMajorFactor> lu;
    lu.compute(normalized);
    if (lu.info() == Success) {
      exponent += M.rows() * scaleExponent;
      return lu.balancedDet(exponent);
    }
    // An aborted factorization met an exactly zero pivot column -- an exactly
    // singular factor -- unless a non-finite entry defeated the pivot search.
    return M.coeffs().allFinite() ? Scalar(0) : Scalar(NumTraits<RealScalar>::quiet_NaN());
  }
};
 
template <typename Factor>
struct kron_factor_ops<Factor, kKronIdentityFactor> {
  // Identity factor, possibly rectangular: the m x n matrix I with I(i,i) = 1
  // for i < min(m, n). I X keeps the leading min(m, n) rows of X, X I^T its
  // leading min(m, n) columns.
  using Scalar = typename Factor::Scalar;
  using TransposedFactor = kron_identity_factor<Scalar, Factor::ColsAtCompileTime, Factor::RowsAtCompileTime>;
  using InverseFactor = Factor;
  static constexpr bool StoresAllEntries = false;
 
  static void prepare(Factor&) {}
  static bool coeffIfStored(const Factor&, Index row, Index col, Scalar& value) {
    if (row != col) return false;
    value = Scalar(1);
    return true;
  }
  template <typename Visitor>
  static void forEachNonZero(const Factor& f, Visitor&& visit) {
    for (Index k = 0; k < numext::mini(f.rows(), f.cols()); ++k) visit(k, k, Scalar(1));
  }
  static Matrix<Index, Dynamic, 1> innerNonZeros(const Factor& f, bool rowMajor) {
    Matrix<Index, Dynamic, 1> counts = Matrix<Index, Dynamic, 1>::Zero(rowMajor ? f.rows() : f.cols());
    counts.head(numext::mini(f.rows(), f.cols())).setOnes();
    return counts;
  }
  static TransposedFactor transposed(const Factor& f) { return TransposedFactor(f.cols(), f.rows()); }
  static const Factor& conjugated(const Factor& f) { return f; }
  static TransposedFactor adjointed(const Factor& f) { return transposed(f); }
  static const Factor& inversed(const Factor& f) { return f; }
  static typename Factor::PlainObject denseFactor(const Factor& f) { return blockOperand(f); }
  static typename Factor::PlainObject::IdentityReturnType blockOperand(const Factor& f) {
    return Factor::PlainObject::Identity(f.rows(), f.cols());
  }
  static bool isSquareIdentity(const Factor& f) { return f.rows() == f.cols(); }
  template <typename Dst, typename Alpha, typename Xpr>
  static void addLeftProduct(Dst& dst, const Alpha& alpha, const Factor& f, const Xpr& X) {
    const Index k = numext::mini(f.rows(), f.cols());
    dst.topRows(k) += alpha * X.topRows(k);
  }
  template <typename Dst, typename Alpha, typename Xpr, typename Work>
  static void addRightProduct(Dst& dst, const Alpha& alpha, const Xpr& X, const Factor& f, Work&) {
    const Index k = numext::mini(f.rows(), f.cols());
    dst.leftCols(k) += alpha * X.leftCols(k);
  }
  static Scalar balancedDet(const Factor&, Index&) { return Scalar(1); }
};
 
template <typename LhsMatrix, typename RhsMatrix>
struct kron_factor_ops<KroneckerOperator<LhsMatrix, RhsMatrix>, kKronKroneckerFactor> {
  // Nested Kronecker factor K = L (x) R: every operation recurses into L and R.
  using Factor = KroneckerOperator<LhsMatrix, RhsMatrix>;
  using Scalar = typename Factor::Scalar;
  using LhsOps = kron_factor_ops<LhsMatrix>;
  using RhsOps = kron_factor_ops<RhsMatrix>;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
  using TransposedFactor = KroneckerOperator<typename LhsOps::TransposedFactor, typename RhsOps::TransposedFactor>;
  using InverseFactor = KroneckerOperator<typename LhsOps::InverseFactor, typename RhsOps::InverseFactor>;
  static constexpr bool StoresAllEntries = LhsOps::StoresAllEntries && RhsOps::StoresAllEntries;
 
  static void prepare(Factor&) {}
  static bool coeffIfStored(const Factor& f, Index row, Index col, Scalar& value) {
    const Index m2 = f.rhs().rows(), n2 = f.rhs().cols();
    Scalar a, b;
    if (!LhsOps::coeffIfStored(f.lhs(), row / m2, col / n2, a) ||
        !RhsOps::coeffIfStored(f.rhs(), row % m2, col % n2, b))
      return false;
    value = a * b;
    return true;
  }
  // Visits (iA m2 + iB, jA n2 + jB) with the loop over L outermost, so a fixed
  // column is visited in increasing row and a fixed row in increasing column,
  // the order the sparse materialization inserts in.
  template <typename Visitor>
  static void forEachNonZero(const Factor& f, Visitor&& visit) {
    using RhsVisitable = kron_factor_visitable<RhsMatrix>;
    const Index m2 = f.rhs().rows(), n2 = f.rhs().cols();
    const auto& R = RhsVisitable::get(f.rhs());
    LhsOps::forEachNonZero(f.lhs(), [&R, &visit, m2, n2](Index iA, Index jA, const Scalar& a) {
      kron_factor_ops<typename RhsVisitable::type>::forEachNonZero(
          R, [&visit, m2, n2, iA, jA, &a](Index iB, Index jB, const Scalar& b) {
            visit(iA * m2 + iB, jA * n2 + jB, a * b);
          });
    });
  }
  static Matrix<Index, Dynamic, 1> innerNonZeros(const Factor& f, bool rowMajor) {
    const Matrix<Index, Dynamic, Dynamic, ColMajor> counts =
        RhsOps::innerNonZeros(f.rhs(), rowMajor) * LhsOps::innerNonZeros(f.lhs(), rowMajor).transpose();
    return counts.reshaped();
  }
  static TransposedFactor transposed(const Factor& f) { return f.transpose(); }
  static Factor conjugated(const Factor& f) { return f.conjugate(); }
  static TransposedFactor adjointed(const Factor& f) { return f.adjoint(); }
  static InverseFactor inversed(const Factor& f) { return f.inverse(); }
  static DenseMatrix denseFactor(const Factor& f) { return DenseMatrix(f); }
  static DenseMatrix blockOperand(const Factor& f) { return denseFactor(f); }
  static bool isSquareIdentity(const Factor&) { return false; }
  template <typename Dst, typename Alpha, typename Xpr>
  static void addLeftProduct(Dst& dst, const Alpha& alpha, const Factor& f, const Xpr& X) {
    f.addProduct(dst, X, alpha);
  }
  // X K^T without forming X^T: with X_j = X(:, j nR : (j+1) nR - 1), j < nL,
  //   X (L (x) R)^T = W L^T,  W(:, j) = vec(X_j R^T),  W (p mR) x nL in work,
  // and W = X_{[p nR x nL]} when R is a square identity.
  template <typename Dst, typename Alpha, typename Xpr, typename Work>
  static void addRightProduct(Dst& dst, const Alpha& alpha, const Xpr& X, const Factor& f, Work& work) {
    const Index p = X.rows(), mL = f.lhs().rows(), nL = f.lhs().cols(), mR = f.rhs().rows(), nR = f.rhs().cols();
    auto dstL = dst.reshaped(p * mR, mL);
    if (RhsOps::isSquareIdentity(f.rhs())) {
      LhsOps::addRightProduct(dstL, alpha, X.reshaped(p * nR, nL), f.lhs(), work);
    } else if (LhsOps::isSquareIdentity(f.lhs())) {
      for (Index j = 0; j < nL; ++j) {
        auto dstj = dst.middleCols(j * mR, mR);
        RhsOps::addRightProduct(dstj, alpha, X.middleCols(j * nR, nR), f.rhs(), work);
      }
    } else {
      Work inner;  // scratch for L and R while work holds W; allocated only if one is nested
      work.setZero(p * mR, nL);
      for (Index j = 0; j < nL; ++j) {
        auto Wj = work.col(j).reshaped(p, mR);
        RhsOps::addRightProduct(Wj, Alpha(1), X.middleCols(j * nR, nR), f.rhs(), inner);
      }
      LhsOps::addRightProduct(dstL, alpha, work, f.lhs(), inner);
    }
  }
  static Scalar balancedDet(const Factor& f, Index& exponent) { return f.balancedDeterminant(exponent); }
};

template <typename Factor, int Kind = kron_factor_kind<Factor>()>
class kron_factor_solver {
 public:
  using Scalar = typename Factor::Scalar;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
 
  explicit kron_factor_solver(const Factor& f) : m_lu(f) {}
  template <typename Xpr>
  DenseMatrix solveLeft(const Xpr& M) const {
    return m_lu.solve(M);
  }
  template <typename Xpr>
  DenseMatrix solveTransposedRight(const Xpr& M) const {
    return m_lu.solve(M.transpose()).transpose();
  }
 
 private:
  PartialPivLU<DenseMatrix> m_lu;
};
 
template <typename Factor>
class kron_factor_solver<Factor, kKronDiagonalFactor> {
 public:
  using Scalar = typename Factor::Scalar;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
 
  explicit kron_factor_solver(const Factor& f) : m_d(f.diagonal()) {}
  template <typename Xpr>
  DenseMatrix solveLeft(const Xpr& M) const {
    return (M.array().colwise() / m_d.array()).matrix();
  }
  template <typename Xpr>
  DenseMatrix solveTransposedRight(const Xpr& M) const {
    return (M.array().rowwise() / m_d.transpose().array()).matrix();
  }
 
 private:
  typename Factor::DiagonalVectorType m_d;
};
 
template <typename Factor>
class kron_factor_solver<Factor, kKronSparseFactor> {
 public:
  using Scalar = typename Factor::Scalar;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
  using ColMajorFactor = SparseMatrix<Scalar, ColMajor, typename Factor::StorageIndex>;
 
  explicit kron_factor_solver(const Factor& f) {
    m_lu.compute(ColMajorFactor(f));
    m_factorized = m_lu.info() == Success;  // else the solves fill with NaN, see kron_factor_ops::nanMatrix
  }
  template <typename Xpr>
  DenseMatrix solveLeft(const Xpr& M) const {
    if (!m_factorized) return kron_factor_ops<Factor>::nanMatrix(M.rows(), M.cols());
    return m_lu.solve(M);
  }
  template <typename Xpr>
  DenseMatrix solveTransposedRight(const Xpr& M) const {
    if (!m_factorized) return kron_factor_ops<Factor>::nanMatrix(M.rows(), M.cols());
    // SparseLU solves into column-major storage only, which a transposed Solve
    // expression would not evaluate into; the right-hand side stays a view.
    const DenseMatrix Xt = m_lu.solve(M.transpose());
    return Xt.transpose();
  }
 
 private:
  bool m_factorized;
  SparseLU<ColMajorFactor> m_lu;
};
 
template <typename Factor>
class kron_factor_solver<Factor, kKronIdentityFactor> {
 public:
  using Scalar = typename Factor::Scalar;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
 
  explicit kron_factor_solver(const Factor&) {}
  template <typename Xpr>
  DenseMatrix solveLeft(const Xpr& M) const {
    return M;
  }
  template <typename Xpr>
  DenseMatrix solveTransposedRight(const Xpr& M) const {
    return M;
  }
};

template <typename LhsMatrix, typename RhsMatrix>
class kron_factor_solver<KroneckerOperator<LhsMatrix, RhsMatrix>, kKronKroneckerFactor> {
 public:
  using Factor = KroneckerOperator<LhsMatrix, RhsMatrix>;
  using Scalar = typename Factor::Scalar;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
 
  explicit kron_factor_solver(const Factor& f)
      : m_solverA(squareFactor(f.lhs())),
        m_solverB(squareFactor(f.rhs())),
        m_n1(f.lhs().cols()),
        m_n2(f.rhs().cols()) {}
  template <typename Xpr>
  DenseMatrix solveLeft(const Xpr& M) const {
    DenseMatrix x(M.rows(), M.cols());
    Factor::solveWith(m_solverA, m_solverB, m_n1, m_n2, M, x);
    return x;
  }
  // M K^{-T} = (K^{-1} M^T)^T.
  template <typename Xpr>
  DenseMatrix solveTransposedRight(const Xpr& M) const {
    return solveLeft(M.transpose()).transpose();
  }
 
 private:
  template <typename F>
  static const F& squareFactor(const F& f) {
    eigen_assert(f.rows() == f.cols() && "KroneckerOperator::solve requires square factors");
    return f;
  }
 
  kron_factor_solver<LhsMatrix> m_solverA;
  kron_factor_solver<RhsMatrix> m_solverB;
  Index m_n1, m_n2;
};

template <typename Factor, int Kind = kron_factor_kind<Factor>()>
struct kron_factor_spectrum {
  using Ops = kron_factor_ops<Factor>;
  using Scalar = typename Factor::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using ComplexScalar = std::complex<RealScalar>;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
  using RealVector = Matrix<RealScalar, Dynamic, 1>;
  using ComplexMatrix = Matrix<ComplexScalar, Dynamic, Dynamic, ColMajor>;
  using ComplexVector = Matrix<ComplexScalar, Dynamic, 1>;
  using Eigenvectors = ComplexMatrix;
  using SingularVectors = DenseMatrix;
 
  static ComplexMatrix complexFactor(const Factor& f) { return Ops::denseFactor(f).template cast<ComplexScalar>(); }
  static ComplexVector eigenvalues(const Factor& f) {
    ComplexEigenSolver<ComplexMatrix> es(complexFactor(f), /*computeEigenvectors=*/false);
    eigen_assert(es.info() == Success);
    return es.eigenvalues();
  }
  static Eigenvectors eigenvectors(const Factor& f) {
    ComplexEigenSolver<ComplexMatrix> es(complexFactor(f));
    eigen_assert(es.info() == Success);
    return es.eigenvectors();
  }
  static RealVector singularValues(const Factor& f) {
    return BDCSVD<DenseMatrix>(Ops::denseFactor(f)).singularValues();
  }
  static SingularVectors matrixU(const Factor& f) {
    return BDCSVD<DenseMatrix, ComputeThinU>(Ops::denseFactor(f)).matrixU();
  }
  static SingularVectors matrixV(const Factor& f) {
    return BDCSVD<DenseMatrix, ComputeThinV>(Ops::denseFactor(f)).matrixV();
  }
};
 
template <typename LhsMatrix, typename RhsMatrix>
struct kron_factor_spectrum<KroneckerOperator<LhsMatrix, RhsMatrix>, kKronKroneckerFactor> {
  using Factor = KroneckerOperator<LhsMatrix, RhsMatrix>;
  using Eigenvectors = KroneckerOperator<typename kron_factor_spectrum<LhsMatrix>::Eigenvectors,
                                         typename kron_factor_spectrum<RhsMatrix>::Eigenvectors>;
  using SingularVectors = KroneckerOperator<typename kron_factor_spectrum<LhsMatrix>::SingularVectors,
                                            typename kron_factor_spectrum<RhsMatrix>::SingularVectors>;
 
  static typename Factor::ComplexVector eigenvalues(const Factor& f) { return f.eigenvalues(); }
  static Eigenvectors eigenvectors(const Factor& f) { return f.eigenvectors(); }
  static typename Factor::RealVector singularValues(const Factor& f) { return f.singularValues(); }
  static SingularVectors matrixU(const Factor& f) { return f.matrixU(); }
  static SingularVectors matrixV(const Factor& f) { return f.matrixV(); }
};
 
}  // namespace internal

template <typename LhsMatrix, typename RhsMatrix>
class KroneckerOperator : public EigenBase<KroneckerOperator<LhsMatrix, RhsMatrix>> {
 public:
  using Scalar = typename LhsMatrix::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using StorageIndex = int;
 
  static_assert(std::is_same<Scalar, typename RhsMatrix::Scalar>::value,
                "KroneckerOperator requires both factors to have the same scalar type");
  static_assert((internal::kron_factor_is_dense_matrix<LhsMatrix>::value ||
                 internal::kron_factor_kind<LhsMatrix>() != internal::kKronDenseFactor) &&
                    (internal::kron_factor_is_dense_matrix<RhsMatrix>::value ||
                     internal::kron_factor_kind<RhsMatrix>() != internal::kKronDenseFactor),
                "KroneckerOperator factors must be plain Matrix, DiagonalMatrix or SparseMatrix types, identity "
                "factors (makeKroneckerOperator stores an Identity() expression as one), KroneckerOperators or "
                "KroneckerSums (owning their storage: views and other expressions would dangle)");
 
 private:
  // Factor-kind dispatch, see kron_factor_ops.
  using LhsOps = internal::kron_factor_ops<LhsMatrix>;
  using RhsOps = internal::kron_factor_ops<RhsMatrix>;
  using LhsSpectrum = internal::kron_factor_spectrum<LhsMatrix>;
  using RhsSpectrum = internal::kron_factor_spectrum<RhsMatrix>;
  template <typename, int>
  friend struct internal::kron_factor_ops;
  template <typename, int>
  friend class internal::kron_factor_solver;
 
 public:
  using ComplexScalar = std::complex<RealScalar>;
  // The vec-trick reshapes below identify a vector of length n1*n2 with an
  // n2 x n1 matrix whose columns are stacked, so every workspace taking part
  // in a reshape is pinned to ColMajor explicitly: the semantics must not
  // change under EIGEN_DEFAULT_TO_ROW_MAJOR.
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
  using RealVector = Matrix<RealScalar, Dynamic, 1>;
  using ComplexMatrix = Matrix<ComplexScalar, Dynamic, Dynamic, ColMajor>;
  using ComplexVector = Matrix<ComplexScalar, Dynamic, 1>;
 
  static constexpr int RowsAtCompileTime =
      internal::size_at_compile_time(LhsMatrix::RowsAtCompileTime, RhsMatrix::RowsAtCompileTime);
  static constexpr int ColsAtCompileTime =
      internal::size_at_compile_time(LhsMatrix::ColsAtCompileTime, RhsMatrix::ColsAtCompileTime);
  static constexpr int MaxRowsAtCompileTime = RowsAtCompileTime;
  static constexpr int MaxColsAtCompileTime = ColsAtCompileTime;
  static constexpr int SizeAtCompileTime = internal::size_at_compile_time(RowsAtCompileTime, ColsAtCompileTime);
  static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime;
  static constexpr bool IsRowMajor = false;
  // Deliberately no IsVectorAtCompileTime: Ref<const KroneckerOperator>'s default
  // StrideType argument reads it, so its absence makes internal::is_ref_compatible
  // SFINAE to false and keeps the iterative solvers on their matrix-free path.

  template <typename LhsDerived, typename RhsDerived>
  KroneckerOperator(const EigenBase<LhsDerived>& a, const EigenBase<RhsDerived>& b)
      : m_A(a.derived()), m_B(b.derived()) {
    eigen_assert(m_A.size() > 0 && m_B.size() > 0 && "KroneckerOperator factors must be non-empty");
    LhsOps::prepare(m_A);
    RhsOps::prepare(m_B);
  }
 
  EIGEN_DEVICE_FUNC Index rows() const { return m_A.rows() * m_B.rows(); }
  EIGEN_DEVICE_FUNC Index cols() const { return m_A.cols() * m_B.cols(); }

  const LhsMatrix& lhs() const { return m_A; }
  const RhsMatrix& rhs() const { return m_B; }

  Scalar coeff(Index row, Index col) const {
    eigen_assert(row >= 0 && row < rows() && col >= 0 && col < cols());
    const Index m2 = m_B.rows(), n2 = m_B.cols();
    Scalar a, b;
    if (!LhsOps::coeffIfStored(m_A, row / m2, col / n2, a) || !RhsOps::coeffIfStored(m_B, row % m2, col % n2, b))
      return Scalar(0);
    return a * b;
  }

  KroneckerOperator<typename LhsOps::TransposedFactor, typename RhsOps::TransposedFactor> transpose() const {
    return {LhsOps::transposed(m_A), RhsOps::transposed(m_B)};
  }

  KroneckerOperator conjugate() const { return {LhsOps::conjugated(m_A), RhsOps::conjugated(m_B)}; }

  KroneckerOperator<typename LhsOps::TransposedFactor, typename RhsOps::TransposedFactor> adjoint() const {
    return {LhsOps::adjointed(m_A), RhsOps::adjointed(m_B)};
  }

  template <typename Rhs>
  Matrix<Scalar, ColsAtCompileTime, Rhs::ColsAtCompileTime> solve(const MatrixBase<Rhs>& b) const {
    EIGEN_STATIC_ASSERT(RowsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic ||
                            int(RowsAtCompileTime) == int(Rhs::RowsAtCompileTime),
                        YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES)
    const Index n1 = m_A.cols(), n2 = m_B.cols();
    eigen_assert(m_A.rows() == n1 && m_B.rows() == n2 && "KroneckerOperator::solve requires square factors");
    eigen_assert(b.rows() == n1 * n2 && "right-hand side has the wrong number of rows");
    const internal::kron_factor_solver<LhsMatrix> solverA(m_A);
    const internal::kron_factor_solver<RhsMatrix> solverB(m_B);
    Matrix<Scalar, ColsAtCompileTime, Rhs::ColsAtCompileTime> x(n1 * n2, b.cols());
    solveWith(solverA, solverB, n1, n2, b.derived(), x);
    return x;
  }

  template <typename Rhs>
  Matrix<Scalar, ColsAtCompileTime, Rhs::ColsAtCompileTime> leastSquaresSolve(const MatrixBase<Rhs>& b) const {
    EIGEN_STATIC_ASSERT(RowsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic ||
                            int(RowsAtCompileTime) == int(Rhs::RowsAtCompileTime),
                        YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES)
    const Index m1 = m_A.rows(), m2 = m_B.rows(), n1 = m_A.cols(), n2 = m_B.cols(), r = b.cols();
    eigen_assert(b.rows() == m1 * m2 && "right-hand side has the wrong number of rows");
    Matrix<Scalar, ColsAtCompileTime, Rhs::ColsAtCompileTime> x(cols(), r);
    // The pivoted QR squares column norms, which over- or underflow on extreme
    // factor magnitudes, so decompose 2^-eA A and 2^-eB B and fold
    // 2^-(eA+eB) back into the solution.
    DenseMatrix An = LhsOps::denseFactor(m_A), Bn = RhsOps::denseFactor(m_B);
    if (!An.allFinite() || !Bn.allFinite()) {
      // The pivoted QR would rank a non-finite column out instead of propagating it.
      x.setConstant(Scalar(NumTraits<RealScalar>::quiet_NaN()));
      return x;
    }
    const int eA = internal::structured_exponent_bound(An), eB = internal::structured_exponent_bound(Bn);
    internal::structured_ldexp_entries_exact(An, -eA);
    internal::structured_ldexp_entries_exact(Bn, -eB);
    const CompleteOrthogonalDecomposition<DenseMatrix> codA(An), codB(Bn);
    typename internal::nested_eval<Rhs, 1>::type actualRhs(b.derived());
    DenseMatrix S, Z, X;
    const Index chunk = internal::kron_rhs_chunk<Scalar>(m1 * m2 + n2 * m1 + n1 * n2);
    for (Index k0 = 0; k0 < r; k0 += chunk) {
      const Index c = numext::mini(chunk, r - k0);
      internal::kron_stack_columns(S, actualRhs.middleCols(k0, c), m2, m1);
      Z = codB.solve(S.reshaped(m2, c * m1));
      X = codA.solve(Z.reshaped(n2 * c, m1).transpose()).transpose();
      internal::structured_ldexp_entries(X, -(eA + eB));
      for (Index k = 0; k < c; ++k) x.col(k0 + k).reshaped(n2, n1) = X.middleRows(k * n2, n2);
    }
    return x;
  }

  Index rank() const {
    BDCSVD<DenseMatrix> svdA(LhsOps::denseFactor(m_A)), svdB(RhsOps::denseFactor(m_B));
    const RealVector sa = svdA.singularValues(), sb = svdB.singularValues();
    // An exactly zero factor zeroes the whole operator (and would make the
    // ratios below 0/0).
    if (sa[0] == RealScalar(0) || sb[0] == RealScalar(0)) return 0;
    const RealScalar tol = relativeRankThreshold();
    const SingularModes modesA(sa), modesB(sb);
    Index r = 0;
    for (Index i = 0; i < sa.size(); ++i)
      for (Index j = 0; j < sb.size(); ++j)
        // Negated ratio comparison so NaN ratios count as non-zero.
        if (modesA.retains(modesB, i, j, tol)) ++r;
    return r;
  }

  KroneckerOperator<typename LhsOps::InverseFactor, typename RhsOps::InverseFactor> inverse() const {
    eigen_assert(m_A.rows() == m_A.cols() && m_B.rows() == m_B.cols() &&
                 "KroneckerOperator::inverse requires square factors");
    return {LhsOps::inversed(m_A), RhsOps::inversed(m_B)};
  }

  Scalar determinant() const {
    eigen_assert(m_A.rows() == m_A.cols() && m_B.rows() == m_B.cols() &&
                 "KroneckerOperator::determinant requires square factors");
    Index exponent = 0;
    const Scalar mant = balancedDeterminant(exponent);
    return internal::structured_ldexp_clamped(mant, exponent);
  }

  ComplexVector eigenvalues() const {
    eigen_assert(m_A.rows() == m_A.cols() && m_B.rows() == m_B.cols() &&
                 "KroneckerOperator::eigenvalues requires square factors");
    // Column-major stacking of the nB x nA outer product puts mu_j(B) lambda_i(A)
    // at index i*nB + j, the Kronecker order.
    return (RhsSpectrum::eigenvalues(m_B) * LhsSpectrum::eigenvalues(m_A).transpose()).reshaped();
  }

  KroneckerOperator<typename LhsSpectrum::Eigenvectors, typename RhsSpectrum::Eigenvectors> eigenvectors() const {
    eigen_assert(m_A.rows() == m_A.cols() && m_B.rows() == m_B.cols() &&
                 "KroneckerOperator::eigenvectors requires square factors");
    return {LhsSpectrum::eigenvectors(m_A), RhsSpectrum::eigenvectors(m_B)};
  }

  RealVector singularValues() const {
    // Column-major stacking of the kB x kA outer product puts sigma_j(B) sigma_i(A)
    // at index i*kB + j, the Kronecker order.
    return (RhsSpectrum::singularValues(m_B) * LhsSpectrum::singularValues(m_A).transpose()).reshaped();
  }

  KroneckerOperator<typename LhsSpectrum::SingularVectors, typename RhsSpectrum::SingularVectors> matrixU() const {
    return {LhsSpectrum::matrixU(m_A), RhsSpectrum::matrixU(m_B)};
  }

  KroneckerOperator<typename LhsSpectrum::SingularVectors, typename RhsSpectrum::SingularVectors> matrixV() const {
    return {LhsSpectrum::matrixV(m_A), RhsSpectrum::matrixV(m_B)};
  }

  template <typename Dest>
  void evalTo(Dest& dst) const {
    evalToImpl(dst, IsSparseDestination<Dest>());
  }

  template <typename Dest>
  void addTo(Dest& dst) const {
    addToImpl(dst, IsSparseDestination<Dest>());
  }

  template <typename Dest>
  void subTo(Dest& dst) const {
    subToImpl(dst, IsSparseDestination<Dest>());
  }
 
 private:
  template <typename Dest>
  using IsSparseDestination = std::is_same<typename internal::traits<Dest>::StorageKind, Sparse>;

  template <typename Dest>
  void evalToImpl(Dest& dst, std::false_type) const {
    const Index m2 = m_B.rows(), n2 = m_B.cols();
    if (!LhsOps::StoresAllEntries) dst.setZero();
    const auto& B = RhsOps::blockOperand(m_B);
    LhsOps::forEachNonZero(
        m_A, [&dst, &B, m2, n2](Index i, Index j, const Scalar& a) { dst.block(i * m2, j * n2, m2, n2) = a * B; });
  }
 
  template <typename Dest>
  void addToImpl(Dest& dst, std::false_type) const {
    const Index m2 = m_B.rows(), n2 = m_B.cols();
    const auto& B = RhsOps::blockOperand(m_B);
    LhsOps::forEachNonZero(
        m_A, [&dst, &B, m2, n2](Index i, Index j, const Scalar& a) { dst.block(i * m2, j * n2, m2, n2) += a * B; });
  }
 
  template <typename Dest>
  void subToImpl(Dest& dst, std::false_type) const {
    const Index m2 = m_B.rows(), n2 = m_B.cols();
    const auto& B = RhsOps::blockOperand(m_B);
    LhsOps::forEachNonZero(
        m_A, [&dst, &B, m2, n2](Index i, Index j, const Scalar& a) { dst.block(i * m2, j * n2, m2, n2) -= a * B; });
  }

  template <typename Dest>
  void evalToImpl(Dest& S, std::true_type) const {
    using LhsVisitable = internal::kron_factor_visitable<LhsMatrix>;
    using RhsVisitable = internal::kron_factor_visitable<RhsMatrix>;
    using VisitedLhsOps = internal::kron_factor_ops<typename LhsVisitable::type>;
    using VisitedRhsOps = internal::kron_factor_ops<typename RhsVisitable::type>;
    const Index m2 = m_B.rows(), n2 = m_B.cols();
    const auto& A = LhsVisitable::get(m_A);
    const auto& B = RhsVisitable::get(m_B);
    S.resize(rows(), cols());
    using IndexVector = Matrix<Index, Dynamic, 1>;
    const IndexVector nnzA = VisitedLhsOps::innerNonZeros(A, Dest::IsRowMajor);
    const IndexVector nnzB = VisitedRhsOps::innerNonZeros(B, Dest::IsRowMajor);
    // Inner vectors kA of A and kB of B meet in inner vector kA * nnzB.size() + kB
    // of the product: the column-major stacking of the count outer product.
    const Matrix<Index, Dynamic, Dynamic, ColMajor> counts = nnzB * nnzA.transpose();
    S.reserve(counts.reshaped());
    VisitedLhsOps::forEachNonZero(A, [&S, &B, m2, n2](Index iA, Index jA, const Scalar& a) {
      VisitedRhsOps::forEachNonZero(B, [&S, m2, n2, iA, jA, &a](Index iB, Index jB, const Scalar& b) {
        S.insert(iA * m2 + iB, jA * n2 + jB) = a * b;
      });
    });
    S.makeCompressed();
  }
 
  template <typename Dest>
  void addToImpl(Dest& dst, std::true_type) const {
    typename Dest::PlainObject product;
    evalTo(product);
    dst += product;
  }
 
  template <typename Dest>
  void subToImpl(Dest& dst, std::true_type) const {
    typename Dest::PlainObject product;
    evalTo(product);
    dst -= product;
  }
 
 public:
  template <typename Rhs>
  Product<KroneckerOperator, Rhs> operator*(const MatrixBase<Rhs>& x) const {
    EIGEN_STATIC_ASSERT(ColsAtCompileTime == Dynamic || Rhs::RowsAtCompileTime == Dynamic ||
                            int(ColsAtCompileTime) == int(Rhs::RowsAtCompileTime),
                        INVALID_MATRIX_PRODUCT)
    eigen_assert(x.rows() == cols() && "invalid product: dimensions do not match");
    return Product<KroneckerOperator, Rhs>(*this, x.derived());
  }

  template <typename Dest, typename Rhs, typename ProductScalar>
  void addProduct(Dest& dst, const Rhs& rhs, const ProductScalar& alpha) const {
    using ProductMatrix = Matrix<ProductScalar, Dynamic, Dynamic, ColMajor>;  // ColMajor: see DenseMatrix
    const Index m1 = m_A.rows(), n1 = m_A.cols(), m2 = m_B.rows(), n2 = m_B.cols(), r = rhs.cols();
    eigen_assert(rhs.rows() == n1 * n2 && "invalid product: dimensions do not match");
    typename internal::nested_eval<Rhs, 1>::type actualRhs(rhs);
    const bool skipA = LhsOps::isSquareIdentity(m_A), skipB = RhsOps::isSquareIdentity(m_B);
    ProductMatrix X, BX, Y, work;
    const Index chunk = internal::kron_rhs_chunk<ProductScalar>(n1 * n2 + m2 * n1 + m2 * m1);
    for (Index k0 = 0; k0 < r; k0 += chunk) {
      const Index c = numext::mini(chunk, r - k0);
      if (c == 1) {
        const auto Xk = actualRhs.col(k0).reshaped(n2, n1);
        auto Yk = dst.col(k0).reshaped(m2, m1);
        if (skipA) {
          RhsOps::addLeftProduct(Yk, alpha, m_B, Xk);
        } else if (skipB) {
          LhsOps::addRightProduct(Yk, alpha, Xk, m_A, work);
        } else {
          BX.setZero(m2, n1);
          RhsOps::addLeftProduct(BX, ProductScalar(1), m_B, Xk);
          LhsOps::addRightProduct(Yk, alpha, BX, m_A, work);
        }
        continue;
      }
      internal::kron_stack_columns(X, actualRhs.middleCols(k0, c), n2, n1);
      if (!skipB) {
        BX.setZero(m2, c * n1);
        RhsOps::addLeftProduct(BX, ProductScalar(1), m_B, X.reshaped(n2, c * n1));
      }
      const auto BXhat = (skipB ? X : BX).reshaped(m2 * c, n1);
      if (!skipA) {
        Y.setZero(m2 * c, m1);
        LhsOps::addRightProduct(Y, ProductScalar(1), BXhat, m_A, work);
      }
      for (Index k = 0; k < c; ++k) {
        auto Yk = dst.col(k0 + k).reshaped(m2, m1);
        if (skipA)
          Yk += alpha * BXhat.middleRows(k * m2, m2);
        else
          Yk += alpha * Y.middleRows(k * m2, m2);
      }
    }
  }
 
 private:
  RealScalar relativeRankThreshold() const {
    return RealScalar(numext::mini(rows(), cols())) * NumTraits<RealScalar>::epsilon();
  }
 
  // Factor-level ratios and frexp decompositions, independent of the other factor.
  struct SingularModes {
    explicit SingularModes(const RealVector& s) : ratios(s.size()), mantissas(s.size()), exponents(s.size()) {
      for (Index i = 0; i < s.size(); ++i) {
        ratios[i] = s[i] / s[0];
        exponents[i] = 0;
        EIGEN_USING_STD(frexp);
        mantissas[i] = frexp(s[i], &exponents[i]);
      }
    }
 
    bool retains(const SingularModes& other, Index i, Index j, RealScalar tol) const {
      // Negation counts NaN ratios as retained, matching SVDBase.
      if (ratios[i] * other.ratios[j] < tol) return false;
      const RealScalar m = mantissas[i] * other.mantissas[j];
      if (!(numext::isfinite)(m)) return true;
      if (m == RealScalar(0)) return false;
      // s_i*s_j = m*2^(e_i+e_j), with m in [0.25,1). Test the
      // smallest-normal clamp without forming an overflowing/underflowing product.
      const int e = exponents[i] + other.exponents[j] - (m < RealScalar(0.5) ? 1 : 0);
      return e >= std::numeric_limits<RealScalar>::min_exponent;
    }
 
    RealVector ratios, mantissas;
    Matrix<int, Dynamic, 1> exponents;
  };

  Scalar balancedDeterminant(Index& exponent) const {
    eigen_assert(m_A.rows() == m_A.cols() && m_B.rows() == m_B.cols() &&
                 "KroneckerOperator::determinant requires square factors");
    const Scalar mant = balancedDetPow(m_A, m_B.cols(), exponent);
    return internal::structured_balance(mant * balancedDetPow(m_B, m_A.cols(), exponent), exponent);
  }

  template <typename SolverA, typename SolverB, typename Rhs, typename Dest>
  static void solveWith(const SolverA& solverA, const SolverB& solverB, Index n1, Index n2, const Rhs& b, Dest& x) {
    typename internal::nested_eval<Rhs, 1>::type actualRhs(b);
    DenseMatrix S, Z, X;
    const Index r = b.cols();
    const Index chunk = internal::kron_rhs_chunk<Scalar>(3 * n1 * n2);
    for (Index k0 = 0; k0 < r; k0 += chunk) {
      const Index c = numext::mini(chunk, r - k0);
      internal::kron_stack_columns(S, actualRhs.middleCols(k0, c), n2, n1);
      Z = solverB.solveLeft(S.reshaped(n2, c * n1));
      X = solverA.solveTransposedRight(Z.reshaped(n2 * c, n1));
      for (Index k = 0; k < c; ++k) x.col(k0 + k).reshaped(n2, n1) = X.middleRows(k * n2, n2);
    }
  }

  template <typename Factor>
  static Scalar balancedDetPow(const Factor& M, Index power, Index& exponent) {
    Index e = 0;
    const Scalar m = internal::kron_factor_ops<Factor>::balancedDet(M, e);
    Scalar r(1);
    Index er = 0;
    for (Index k = 0; k < power; ++k) r = internal::structured_balance(r * m, er);
    exponent += power * e + er;
    return r;
  }
 
  LhsMatrix m_A;
  RhsMatrix m_B;
};

template <typename LhsDerived, typename RhsDerived>
KroneckerOperator<typename internal::kron_factor_storage<LhsDerived>::type,
                  typename internal::kron_factor_storage<RhsDerived>::type>
makeKroneckerOperator(const EigenBase<LhsDerived>& a, const EigenBase<RhsDerived>& b) {
  return {a.derived(), b.derived()};
}

template <typename D1, typename D2, typename D3, typename... Rest>
auto makeKroneckerOperator(const EigenBase<D1>& a, const EigenBase<D2>& b, const EigenBase<D3>& c,
                           const Rest&... rest) {
  return makeKroneckerOperator(a, makeKroneckerOperator(b, c, rest...));
}
 
namespace internal {
 
template <typename LhsMatrix, typename RhsMatrix, typename Rhs, int ProductTag>
struct generic_product_impl<KroneckerOperator<LhsMatrix, RhsMatrix>, Rhs, StructuredShape, DenseShape, ProductTag>
    : structured_product_impl<KroneckerOperator<LhsMatrix, RhsMatrix>, Rhs> {};
 
}  // namespace internal
 
}  // namespace Eigen
 
#endif  // EIGEN_STRUCTURED_KRONECKER_OPERATOR_H