KroneckerSum.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] R. H. Bartels and G. W. Stewart, "Solution of the matrix equation
//      AX + XB = C", Communications of the ACM 15 (1972), 820-826. The Schur
//      form solve of the Sylvester equation B X + X A^T = mat(b) that
//      BartelsStewart generalizes to Kronecker sums of any number of factors.
//  [2] R. E. Lynch, J. R. Rice and D. H. Thomas, "Direct solution of partial
//      difference equations by tensor product methods", Numerische Mathematik 6
//      (1964), 185-199. The fast diagonalization method: with Hermitian factors
//      the Schur forms are diagonal and the solve is a division by the
//      eigenvalue sums.
 
#ifndef EIGEN_STRUCTURED_KRONECKER_SUM_H
#define EIGEN_STRUCTURED_KRONECKER_SUM_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
template <typename KroneckerSumType>
class BartelsStewart;
 
namespace internal {
 
template <typename LhsMatrix, typename RhsMatrix>
struct traits<KroneckerSum<LhsMatrix, RhsMatrix>> {
  using Scalar = typename LhsMatrix::Scalar;
  using StorageKind = Dense;
  using XprKind = MatrixXpr;
  using StorageIndex = int;
  static constexpr int RowsAtCompileTime =
      size_at_compile_time(traits<LhsMatrix>::RowsAtCompileTime, traits<RhsMatrix>::RowsAtCompileTime);
  static constexpr int ColsAtCompileTime = RowsAtCompileTime;
  static constexpr int MaxRowsAtCompileTime = RowsAtCompileTime;
  static constexpr int MaxColsAtCompileTime = RowsAtCompileTime;
  // No NestByRefBit, for the reason given at traits<KroneckerOperator>.
  static constexpr unsigned int Flags = 0;
};
 
template <typename LhsMatrix, typename RhsMatrix>
struct evaluator_traits<KroneckerSum<LhsMatrix, RhsMatrix>> {
  using Kind = IndexBased;
  using Shape = StructuredShape;
};
 
template <typename KroneckerSumType>
struct traits<BartelsStewart<KroneckerSumType>> : traits<Matrix<typename KroneckerSumType::Scalar, Dynamic, Dynamic>> {
  using XprKind = MatrixXpr;
  using StorageKind = SolverStorage;
  using StorageIndex = int;
  using BaseTraits = traits<Matrix<typename KroneckerSumType::Scalar, Dynamic, Dynamic>>;
  static constexpr unsigned int Flags = BaseTraits::Flags & RowMajorBit;
  static constexpr int CoeffReadCost = Dynamic;
};

template <typename LhsMatrix, typename RhsMatrix>
struct kron_factor_ops<KroneckerSum<LhsMatrix, RhsMatrix>, kKronSumFactor> {
  using Factor = KroneckerSum<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 SparseFactor = SparseMatrix<Scalar>;
  using TransposedFactor = KroneckerSum<typename LhsOps::TransposedFactor, typename RhsOps::TransposedFactor>;
  using InverseFactor = DenseMatrix;
  static constexpr bool StoresAllEntries = false;
 
  static void prepare(Factor&) {}
  static bool coeffIfStored(const Factor& f, Index row, Index col, Scalar& value) {
    const Index n2 = f.rhs().rows();
    Scalar a(0), b(0);
    const bool storedA = row % n2 == col % n2 && LhsOps::coeffIfStored(f.lhs(), row / n2, col / n2, a);
    const bool storedB = row / n2 == col / n2 && RhsOps::coeffIfStored(f.rhs(), row % n2, col % n2, b);
    if (!storedA && !storedB) return false;
    value = a + b;
    return true;
  }
  template <typename Visitor>
  static void forEachNonZero(const Factor& f, Visitor&& visit) {
    kron_factor_ops<SparseFactor>::forEachNonZero(kron_factor_visitable<Factor>::get(f), visit);
  }
  static Matrix<Index, Dynamic, 1> innerNonZeros(const Factor& f, bool rowMajor) {
    return kron_factor_ops<SparseFactor>::innerNonZeros(kron_factor_visitable<Factor>::get(f), rowMajor);
  }
  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.solve(DenseMatrix::Identity(f.rows(), f.cols())); }
  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);
  }
  // With X_j = X(:, j n2 : (j+1) n2 - 1), j < n1, the p x n2 slices of X,
  //   X (L (+) R)^T = X (L (x) I)^T + X (I (x) R)^T = X_{[p n2 x n1]} L^T + [X_j R^T]_j.
  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(), n1 = f.lhs().rows(), n2 = f.rhs().rows();
    auto dstL = dst.reshaped(p * n2, n1);
    LhsOps::addRightProduct(dstL, alpha, X.reshaped(p * n2, n1), f.lhs(), work);
    for (Index j = 0; j < n1; ++j) {
      auto dstj = dst.middleCols(j * n2, n2);
      RhsOps::addRightProduct(dstj, alpha, X.middleCols(j * n2, n2), f.rhs(), work);
    }
  }
  // det = prod_{i,j} (lambda_i + mu_j) would cost only the factor Schur forms,
  // but every delta lambda_i recurs in n_R factors of the product, and a shift
  // split across L and R makes |lambda_i + mu_j| << |lambda_i| + |mu_j|. Against
  // a 256-bit reference that product lost 1 to 6 digits to the LU of the
  // materialized sum.
  static Scalar balancedDet(const Factor& f, Index& exponent) {
    return kron_factor_ops<DenseMatrix>::balancedDet(denseFactor(f), exponent);
  }
};
 
template <typename LhsMatrix, typename RhsMatrix>
struct kron_factor_visitable<KroneckerSum<LhsMatrix, RhsMatrix>, kKronSumFactor> {
  using Factor = KroneckerSum<LhsMatrix, RhsMatrix>;
  using type = SparseMatrix<typename Factor::Scalar>;
  static type get(const Factor& f) {
    type S;
    S = f;
    return S;
  }
};

template <typename LhsMatrix, typename RhsMatrix>
struct kron_factor_spectrum<KroneckerSum<LhsMatrix, RhsMatrix>, kKronSumFactor>
    : kron_factor_spectrum<KroneckerSum<LhsMatrix, RhsMatrix>, kKronDenseFactor> {
  using Factor = KroneckerSum<LhsMatrix, RhsMatrix>;
  using Eigenvectors = KroneckerOperator<typename kron_factor_spectrum<LhsMatrix>::Eigenvectors,
                                         typename kron_factor_spectrum<RhsMatrix>::Eigenvectors>;
 
  static typename Factor::ComplexVector eigenvalues(const Factor& f) { return f.eigenvalues(); }
  static Eigenvectors eigenvectors(const Factor& f) { return f.eigenvectors(); }
};
 
template <typename LhsMatrix, typename RhsMatrix>
class kron_factor_solver<KroneckerSum<LhsMatrix, RhsMatrix>, kKronSumFactor> {
 public:
  using Factor = KroneckerSum<LhsMatrix, RhsMatrix>;
  using Scalar = typename Factor::Scalar;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
 
  explicit kron_factor_solver(const Factor& f) : m_solver(f) {}
  template <typename Xpr>
  DenseMatrix solveLeft(const Xpr& M) const {
    return m_solver.solve(M);
  }
  // X = M S^{-T} solves S X^T = M^T, written through a transposed view of X.
  template <typename Xpr>
  DenseMatrix solveTransposedRight(const Xpr& M) const {
    DenseMatrix X(M.rows(), M.cols());
    X.transpose() = m_solver.solve(M.transpose());
    return X;
  }
 
 private:
  BartelsStewart<Factor> m_solver;
};
 
}  // namespace internal

template <typename LhsMatrix, typename RhsMatrix>
class KroneckerSum : public EigenBase<KroneckerSum<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,
                "KroneckerSum 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),
                "KroneckerSum factors must be plain Matrix, DiagonalMatrix or SparseMatrix types, identity factors "
                "(makeKroneckerSum stores an Identity() expression as one), KroneckerOperators or KroneckerSums "
                "(owning their storage)");
 
 private:
  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>;
 
 public:
  using ComplexScalar = std::complex<RealScalar>;
  using DenseMatrix = Matrix<Scalar, 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 = RowsAtCompileTime;
  static constexpr int MaxRowsAtCompileTime = RowsAtCompileTime;
  static constexpr int MaxColsAtCompileTime = RowsAtCompileTime;
  static constexpr int SizeAtCompileTime = internal::size_at_compile_time(RowsAtCompileTime, ColsAtCompileTime);
  static constexpr int MaxSizeAtCompileTime = SizeAtCompileTime;
  static constexpr bool IsRowMajor = false;
  // Deliberately no IsVectorAtCompileTime, for the reason given in KroneckerOperator.

  template <typename LhsDerived, typename RhsDerived>
  KroneckerSum(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 && "KroneckerSum factors must be non-empty");
    eigen_assert(m_A.rows() == m_A.cols() && m_B.rows() == m_B.cols() && "KroneckerSum factors must be square");
    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 rows(); }

  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());
    Scalar value;
    return internal::kron_factor_ops<KroneckerSum>::coeffIfStored(*this, row, col, value) ? value : Scalar(0);
  }

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

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

  KroneckerSum<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 {
    const BartelsStewart<KroneckerSum> solver(*this);
    return solver.solve(b);
  }

  ComplexVector eigenvalues() const {
    const ComplexVector lambda = LhsSpectrum::eigenvalues(m_A), mu = RhsSpectrum::eigenvalues(m_B);
    return (mu.replicate(fix<1>, lambda.size()) + lambda.transpose().replicate(mu.size(), fix<1>)).reshaped();
  }

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

  template <typename Rhs>
  Product<KroneckerSum, 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<KroneckerSum, 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>;
    const Index n1 = m_A.rows(), n2 = m_B.rows(), r = rhs.cols();
    eigen_assert(rhs.rows() == n1 * n2 && "invalid product: dimensions do not match");
    typename internal::nested_eval<Rhs, 1>::type actualRhs(rhs);
    ProductMatrix X, Y, work;
    const Index chunk = internal::kron_rhs_chunk<ProductScalar>(2 * n1 * n2);
    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(n2, n1);
        RhsOps::addLeftProduct(Yk, alpha, m_B, Xk);
        LhsOps::addRightProduct(Yk, alpha, Xk, m_A, work);
        continue;
      }
      internal::kron_stack_columns(X, actualRhs.middleCols(k0, c), n2, n1);
      Y.setZero(n2 * c, n1);
      auto Yflat = Y.reshaped(n2, c * n1);
      RhsOps::addLeftProduct(Yflat, ProductScalar(1), m_B, X.reshaped(n2, c * n1));
      LhsOps::addRightProduct(Y, ProductScalar(1), X, m_A, work);
      for (Index k = 0; k < c; ++k) dst.col(k0 + k).reshaped(n2, n1) += alpha * Y.middleRows(k * n2, n2);
    }
  }

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

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

  template <typename Dest>
  void subTo(Dest& dst) const {
    addToImpl(dst, Scalar(-1), 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 {
    dst.setZero();
    addToImpl(dst, Scalar(1), std::false_type());
  }

  template <typename Dest>
  void addToImpl(Dest& dst, const Scalar& s, std::false_type) const {
    const Index n1 = m_A.rows(), n2 = m_B.rows();
    LhsOps::forEachNonZero(m_A, [&dst, &s, n2](Index i, Index j, const Scalar& a) {
      dst.block(i * n2, j * n2, n2, n2).diagonal().array() += s * a;
    });
    const auto& B = RhsOps::blockOperand(m_B);
    for (Index k = 0; k < n1; ++k) dst.block(k * n2, k * n2, n2, n2) += s * B;
  }

  template <typename Dest>
  void evalToImpl(Dest& S, std::true_type) const {
    using IndexVector = Matrix<Index, Dynamic, 1>;
    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 n1 = m_A.rows(), n2 = m_B.rows();
    const auto& A = LhsVisitable::get(m_A);
    const auto& B = RhsVisitable::get(m_B);
    const IndexVector nnzA = VisitedLhsOps::innerNonZeros(A, Dest::IsRowMajor);
    const IndexVector nnzB = VisitedRhsOps::innerNonZeros(B, Dest::IsRowMajor);
    Dest SA(rows(), cols()), SB(rows(), cols());
    SA.reserve(IndexVector(nnzA.transpose().replicate(n2, fix<1>).reshaped()));
    VisitedLhsOps::forEachNonZero(A, [&SA, n2](Index i, Index j, const Scalar& a) {
      for (Index k = 0; k < n2; ++k) SA.insert(i * n2 + k, j * n2 + k) = a;
    });
    SB.reserve(IndexVector(nnzB.replicate(n1, fix<1>)));
    for (Index k = 0; k < n1; ++k)
      VisitedRhsOps::forEachNonZero(
          B, [&SB, n2, k](Index i, Index j, const Scalar& b) { SB.insert(k * n2 + i, k * n2 + j) = b; });
    S = SA + SB;
  }
 
  template <typename Dest>
  void addToImpl(Dest& dst, const Scalar& s, std::true_type) const {
    typename Dest::PlainObject sum;
    evalTo(sum);
    if (s == Scalar(1))
      dst += sum;
    else
      dst -= sum;
  }
 
  LhsMatrix m_A;
  RhsMatrix m_B;
};

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

template <typename D1, typename D2, typename D3, typename... Rest>
auto makeKroneckerSum(const EigenBase<D1>& a, const EigenBase<D2>& b, const EigenBase<D3>& c, const Rest&... rest) {
  return makeKroneckerSum(a, makeKroneckerSum(b, c, rest...));
}

template <typename KroneckerSumType>
class BartelsStewart : public SolverBase<BartelsStewart<KroneckerSumType>> {
 public:
  using Base = SolverBase<BartelsStewart>;
  friend class SolverBase<BartelsStewart>;
  EIGEN_GENERIC_PUBLIC_INTERFACE(BartelsStewart)
  EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar)
  using MatrixType = KroneckerSumType;
  using ComplexScalar = std::complex<RealScalar>;
  using DenseMatrix = Matrix<Scalar, Dynamic, Dynamic, ColMajor>;
  using ComplexMatrix = Matrix<ComplexScalar, Dynamic, Dynamic, ColMajor>;
  using RealVector = Matrix<RealScalar, Dynamic, 1>;

  BartelsStewart() = default;

  explicit BartelsStewart(const KroneckerSumType& op) { compute(op); }

  BartelsStewart& compute(const KroneckerSumType& op) {
    std::vector<DenseMatrix> leaves;
    collectLeaves(op, leaves);
    const std::size_t d = leaves.size();
    m_sizes.resize(d);
    m_inner.resize(d);
    m_unitary.clear();
    m_triangular.clear();
    m_triangularAdjoint.clear();
    m_basis.clear();
    m_quasiTriangular.clear();
    m_info = Success;
    bool hermitian = true;
    for (std::size_t k = 0; k < d; ++k) {
      m_sizes[k] = leaves[k].rows();
      if (!leaves[k].allFinite()) m_info = InvalidInput;
      hermitian = hermitian && leaves[k] == leaves[k].adjoint();
    }
    if (hermitian)
      m_path = Path::Diagonal;
    else if (!NumTraits<Scalar>::IsComplex && d <= std::size_t(kMaxRealFactors))
      m_path = Path::RealSchur;
    else
      m_path = Path::ComplexSchur;
    m_size = 1;
    for (std::size_t k = d; k-- > 0;) {
      m_inner[k] = m_size;
      m_size *= m_sizes[k];
    }
    m_isInitialized = true;
    if (m_info != Success) return *this;
    switch (m_path) {
      case Path::Diagonal: {
        m_spectrum = RealVector::Zero(1);
        SelfAdjointEigenSolver<DenseMatrix> es;
        for (std::size_t k = 0; k < d; ++k) {
          es.compute(leaves[k]);
          if (es.info() != Success) m_info = NoConvergence;
          m_basis.push_back(es.eigenvectors());
          // Kronecker order: the new factor's index runs fastest.
          const RealVector& lambda = es.eigenvalues();
          const RealVector previous = m_spectrum;
          m_spectrum =
              (lambda.replicate(fix<1>, previous.size()) + previous.transpose().replicate(lambda.size(), fix<1>))
                  .reshaped();
        }
        break;
      }
      case Path::RealSchur:
        computeRealSchur(leaves, internal::bool_constant<!NumTraits<Scalar>::IsComplex>());
        break;
      case Path::ComplexSchur: {
        ComplexSchur<ComplexMatrix> schur;
        for (std::size_t k = 0; k < d; ++k) {
          schur.compute(leaves[k].template cast<ComplexScalar>());
          if (schur.info() != Success) m_info = NoConvergence;
          m_unitary.push_back(schur.matrixU());
          m_triangular.push_back(schur.matrixT());
          m_triangularAdjoint.emplace_back(m_triangular.back().adjoint());
        }
        break;
      }
    }
    return *this;
  }
 
  Index rows() const noexcept { return m_size; }
  Index cols() const noexcept { return m_size; }

  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "BartelsStewart is not initialized.");
    return m_info;
  }

  bool isHermitian() const {
    eigen_assert(m_isInitialized && "BartelsStewart is not initialized.");
    return m_path == Path::Diagonal;
  }
 
#ifdef EIGEN_PARSED_BY_DOXYGEN
  template <typename Rhs>
  inline const Solve<BartelsStewart, Rhs> solve(const MatrixBase<Rhs>& b) const;
#endif
 
#ifndef EIGEN_PARSED_BY_DOXYGEN
  template <typename RhsType, typename DstType>
  void _solve_impl(const RhsType& rhs, DstType& dst) const {
    solveImpl(rhs, dst, /*adjoint=*/false);
  }
 
  // M^T = conj(M^H): M^{-T} b = conj(M^{-H} conj(b)).
  template <bool Conjugate, typename RhsType, typename DstType>
  void _solve_impl_transposed(const RhsType& rhs, DstType& dst) const {
    constexpr bool ConjugateRhs = !Conjugate && NumTraits<Scalar>::IsComplex;
    solveImpl(rhs.template conjugateIf<ConjugateRhs>(), dst, /*adjoint=*/true);
    if (ConjugateRhs) dst = dst.conjugate();
  }
#endif
 
 private:
  template <typename RhsType, typename DstType>
  void solveImpl(const RhsType& rhs, DstType& dst, bool adjoint) const {
    if (m_info != Success) {
      // No usable decompositions; the nested solves never see info().
      dst.setConstant(Scalar(NumTraits<RealScalar>::quiet_NaN()));
      return;
    }
    switch (m_path) {
      case Path::Diagonal: {
        DenseMatrix W = rhs;
        applyBasis(W, m_basis, /*adjoint=*/true);
        W.array().colwise() /= m_spectrum.array();
        applyBasis(W, m_basis, /*adjoint=*/false);
        dst = W;
        break;
      }
      case Path::RealSchur: {
        DenseMatrix W = rhs;
        solveRealSchur(W, adjoint, internal::bool_constant<!NumTraits<Scalar>::IsComplex>());
        dst = W;
        break;
      }
      case Path::ComplexSchur: {
        ComplexMatrix W = rhs.template cast<ComplexScalar>();
        applyBasis(W, m_unitary, /*adjoint=*/true);
        for (Index j = 0; j < W.cols(); ++j) {
          if (adjoint)
            adjointTriangularSolve(0, ComplexScalar(0), W.col(j).data());
          else
            triangularSolve(0, ComplexScalar(0), W.col(j).data());
        }
        applyBasis(W, m_unitary, /*adjoint=*/false);
        dst = internal::structured_scalar_part_impl<Scalar>::run(W);
        break;
      }
    }
  }
 
  template <typename Factor>
  static void collectLeaves(const Factor& f, std::vector<DenseMatrix>& leaves) {
    collectLeaves(f, leaves, internal::kron_factor_is_kronecker_sum<Factor>());
  }
  template <typename Factor>
  static void collectLeaves(const Factor& f, std::vector<DenseMatrix>& leaves, std::true_type) {
    collectLeaves(f.lhs(), leaves);
    collectLeaves(f.rhs(), leaves);
  }
  template <typename Factor>
  static void collectLeaves(const Factor& f, std::vector<DenseMatrix>& leaves, std::false_type) {
    leaves.push_back(DenseMatrix(internal::kron_factor_ops<Factor>::denseFactor(f)));
  }

  template <typename Work, typename Basis>
  void applyBasis(Work& W, const std::vector<Basis>& Q, bool adjoint) const {
    using WorkMatrix = Matrix<typename Work::Scalar, Dynamic, Dynamic, ColMajor>;
    WorkMatrix T;
    for (std::size_t k = 0; k < Q.size(); ++k) {
      const Index n = m_sizes[k], s = m_inner[k], blocks = W.size() / (n * s);
      if (s == 1) {
        Map<WorkMatrix> U(W.data(), n, blocks);
        if (adjoint)
          T.noalias() = Q[k].adjoint() * U;
        else
          T.noalias() = Q[k] * U;
        U = T;
        continue;
      }
      for (Index b = 0; b < blocks; ++b) {
        Map<WorkMatrix> V(W.data() + b * n * s, s, n);
        if (adjoint)
          T.noalias() = V * Q[k].conjugate();
        else
          T.noalias() = V * Q[k].transpose();
        V = T;
      }
    }
  }

  void triangularSolve(std::size_t k, const ComplexScalar& sigma, ComplexScalar* y) const {
    const Index n = m_sizes[k], s = m_inner[k];
    const TriangularMatrix& T = m_triangular[k];
    if (k + 1 == m_sizes.size()) {
      // Last factor, s = 1: back substitution on the shifted T_d, by rows.
      Map<Matrix<ComplexScalar, 1, Dynamic>> Y(y, n);
      for (Index i = n - 1; i >= 0; --i) {
        const Index tail = n - 1 - i;
        Y(i) = (Y(i) - Y.tail(tail).cwiseProduct(T.row(i).tail(tail)).sum()) / (sigma + T(i, i));
      }
      return;
    }
    Map<ComplexMatrix> Y(y, s, n);
    for (Index i = n - 1; i >= 0; --i) {
      const Index tail = n - 1 - i;
      if (tail > 0) Y.col(i).noalias() -= Y.rightCols(tail) * T.row(i).tail(tail).transpose();
      triangularSolve(k + 1, sigma + T(i, i), Y.col(i).data());
    }
  }

  void adjointTriangularSolve(std::size_t k, const ComplexScalar& sigma, ComplexScalar* y) const {
    const Index n = m_sizes[k], s = m_inner[k];
    const TriangularMatrix& L = m_triangularAdjoint[k];
    if (k + 1 == m_sizes.size()) {
      // Last factor, s = 1: forward substitution on the shifted L_d, by rows.
      Map<Matrix<ComplexScalar, 1, Dynamic>> Y(y, n);
      for (Index i = 0; i < n; ++i) Y(i) = (Y(i) - Y.head(i).cwiseProduct(L.row(i).head(i)).sum()) / (sigma + L(i, i));
      return;
    }
    Map<ComplexMatrix> Y(y, s, n);
    for (Index i = 0; i < n; ++i) {
      if (i > 0) Y.col(i).noalias() -= Y.leftCols(i) * L.row(i).head(i).transpose();
      adjointTriangularSolve(k + 1, sigma + L(i, i), Y.col(i).data());
    }
  }
 
  // The real path's coupled systems are at most 2^d wide, see the class documentation.
  static constexpr int kMaxRealFactors = 4;
  static constexpr int kMaxCoupling = 1 << kMaxRealFactors;
  using CouplingMatrix = Matrix<RealScalar, Dynamic, Dynamic, ColMajor, kMaxCoupling, kMaxCoupling>;
  using CouplingVector = Matrix<RealScalar, Dynamic, 1, ColMajor, kMaxCoupling, 1>;
  using RealMatrix = Matrix<RealScalar, Dynamic, Dynamic, ColMajor>;
 
  void computeRealSchur(const std::vector<DenseMatrix>& leaves, std::true_type) {
    RealSchur<DenseMatrix> schur;
    for (const DenseMatrix& leaf : leaves) {
      schur.compute(leaf);
      if (schur.info() != Success) m_info = NoConvergence;
      m_basis.push_back(schur.matrixU());
      m_quasiTriangular.push_back(schur.matrixT());
    }
  }
  void computeRealSchur(const std::vector<DenseMatrix>&, std::false_type) {}
 
  void solveRealSchur(DenseMatrix& W, bool adjoint, std::true_type) const {
    applyBasis(W, m_basis, /*adjoint=*/true);
    std::vector<RealMatrix> work(m_sizes.size());
    const CouplingMatrix outer = CouplingMatrix::Zero(1, 1);
    for (Index j = 0; j < W.cols(); ++j)
      quasiTriangularSolve(0, outer, Map<RealMatrix>(W.col(j).data(), W.rows(), 1), work, adjoint);
    applyBasis(W, m_basis, /*adjoint=*/false);
  }
  void solveRealSchur(DenseMatrix&, bool, std::false_type) const {}

  void quasiTriangularSolve(std::size_t k, const CouplingMatrix& Sigma, Map<RealMatrix> X,
                            std::vector<RealMatrix>& work, bool adjoint) const {
    const Index n = m_sizes[k], s = m_inner[k], m = Sigma.rows();
    const QuasiTriangularMatrix& T = m_quasiTriangular[k];
    const bool last = k + 1 == m_sizes.size();
    for (Index done = 0; done < n;) {
      // Block I = [i, i + p): the next one up from the bottom, or down from the top when adjoint.
      Index i, p;
      if (adjoint) {
        i = done;
        p = i + 1 < n && T(i + 1, i) != RealScalar(0) ? 2 : 1;
      } else {
        const Index end = n - done;
        p = end > 1 && T(end - 1, end - 2) != RealScalar(0) ? 2 : 1;
        i = end - p;
      }
      const Index tail = n - i - p;
      // T(I,I), T(J,I) for J < I, and T(I,J) for J > I.
      const auto TII = T.block(i, i, p, p);
      const auto colHead = T.block(0, i, i, p);
      const auto rowTail = T.block(i, i + p, p, tail);
      const CouplingMatrix coupled = adjoint ? kroneckerSum(Sigma, TII.transpose()) : kroneckerSum(Sigma, TII);
      if (last) {
        // s = 1: row i of X belongs to index i of the last factor.
        auto XI = X.middleRows(i, p);
        if (adjoint) {
          if (i > 0) XI.noalias() -= colHead.transpose().lazyProduct(X.topRows(i));
        } else if (tail > 0) {
          XI.noalias() -= rowTail.lazyProduct(X.bottomRows(tail));
        }
        solveCoupled(coupled, XI);
      } else {
        for (Index t = 0; t < m; ++t) {
          auto Xt = X.col(t).reshaped(s, n);
          auto XtI = Xt.middleCols(i, p);
          if (adjoint) {
            if (i > 0) XtI.noalias() -= Xt.leftCols(i) * colHead;
          } else if (tail > 0) {
            XtI.noalias() -= Xt.rightCols(tail) * rowTail.transpose();
          }
        }
        if (m == 1) {
          quasiTriangularSolve(k + 1, coupled, Map<RealMatrix>(X.col(0).data() + i * s, s, p), work, adjoint);
        } else {
          // Columns (t, a) of Z, a fastest: column a of block I for coupling index t.
          RealMatrix& Z = work[k];
          Z.resize(s, m * p);
          Z.reshaped(p * s, m) = X.middleRows(i * s, p * s);
          quasiTriangularSolve(k + 1, coupled, Map<RealMatrix>(Z.data(), s, m * p), work, adjoint);
          X.middleRows(i * s, p * s) = Z.reshaped(p * s, m);
        }
      }
      done += p;
    }
  }

  template <typename Block>
  static CouplingMatrix kroneckerSum(const CouplingMatrix& Sigma, const Block& S) {
    const Index m = Sigma.rows(), p = S.rows();
    eigen_internal_assert(m * p <= kMaxCoupling);
    CouplingMatrix R = CouplingMatrix::Zero(m * p, m * p);
    for (Index t = 0; t < m; ++t) {
      for (Index u = 0; u < m; ++u) R.block(t * p, u * p, p, p).diagonal().setConstant(Sigma(t, u));
      R.block(t * p, t * p, p, p) += S;
    }
    return R;
  }

  template <typename Block>
  static void solveCoupled(const CouplingMatrix& coupled, Block& XI) {
    if (coupled.rows() == 1) {
      XI(0, 0) /= coupled(0, 0);
      return;
    }
    CouplingVector z = XI.reshaped();
    z = coupled.partialPivLu().solve(z);
    XI.reshaped() = z;
  }
 
  std::vector<Index> m_sizes, m_inner;
  Index m_size = 0;
  std::vector<DenseMatrix> m_basis;
  RealVector m_spectrum;
  // Row-major: the back substitution reads T_k by rows, and the forward
  // substitution of the adjoint solve reads T_k^H, kept as its own copy, the
  // same way.
  using TriangularMatrix = Matrix<ComplexScalar, Dynamic, Dynamic, RowMajor>;
  using QuasiTriangularMatrix = Matrix<RealScalar, Dynamic, Dynamic, RowMajor>;
  std::vector<ComplexMatrix> m_unitary;
  std::vector<TriangularMatrix> m_triangular, m_triangularAdjoint;
  std::vector<QuasiTriangularMatrix> m_quasiTriangular;
  enum class Path { Diagonal, RealSchur, ComplexSchur };
  Path m_path = Path::Diagonal;
  bool m_isInitialized = false;
  ComputationInfo m_info = InvalidInput;
};
 
namespace internal {
 
template <typename LhsMatrix, typename RhsMatrix, typename Rhs, int ProductTag>
struct generic_product_impl<KroneckerSum<LhsMatrix, RhsMatrix>, Rhs, StructuredShape, DenseShape, ProductTag>
    : structured_product_impl<KroneckerSum<LhsMatrix, RhsMatrix>, Rhs> {};
 
}  // namespace internal
 
}  // namespace Eigen
 
#endif  // EIGEN_STRUCTURED_KRONECKER_SUM_H