TridiagonalEigenSolver.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2026 Rasmus Munk Larsen <rmlarsen@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_TRIDIAGONAL_EIGENSOLVER_H
#define EIGEN_TRIDIAGONAL_EIGENSOLVER_H
 
#include "./TridiagonalBisection.h"
#include "./TridiagonalInverseIteration.h"
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {

template <typename Scalar_>
class TridiagonalEigenSolver {
 public:
  using Scalar = Scalar_;
  using RealScalar = Scalar;
  static_assert(NumTraits<Scalar>::IsComplex == 0 && NumTraits<Scalar>::IsInteger == 0,
                "TridiagonalEigenSolver requires a real floating-point scalar type; for the eigenvalues of a complex "
                "Hermitian tridiagonal matrix, pass subdiag.cwiseAbs() as the real sub-diagonal (see the class "
                "documentation).");

  using VectorType = Matrix<Scalar, Dynamic, 1>;
  using MatrixType = Matrix<Scalar, Dynamic, Dynamic>;

  TridiagonalEigenSolver() = default;

  explicit TridiagonalEigenSolver(Index size)
      : m_eivalues(size), m_eivec(size, size), m_diag(size), m_subdiag(size > 1 ? size - 1 : 0) {}

  template <typename DiagType, typename SubdiagType>
  TridiagonalEigenSolver(const MatrixBase<DiagType>& diag, const MatrixBase<SubdiagType>& subdiag,
                         int options = ComputeEigenvectors, const EigenvalueRange& range = EigenvalueRange::all())
      : TridiagonalEigenSolver() {
    compute(diag, subdiag, options, range);
  }

  template <typename DiagType, typename SubdiagType>
  TridiagonalEigenSolver& compute(const MatrixBase<DiagType>& diag, const MatrixBase<SubdiagType>& subdiag,
                                  int options = ComputeEigenvectors,
                                  const EigenvalueRange& range = EigenvalueRange::all()) {
    eigen_assert((options & ~EigVecMask) == 0 && (options & EigVecMask) != EigVecMask && "invalid option parameter");
    computeEigenvalues(diag, subdiag, range);
    if (m_info == Success && (options & ComputeEigenvectors) == ComputeEigenvectors) computeEigenvectors();
    return *this;
  }

  template <typename DiagType, typename SubdiagType>
  TridiagonalEigenSolver& computeEigenvalues(const MatrixBase<DiagType>& diag, const MatrixBase<SubdiagType>& subdiag,
                                             const EigenvalueRange& range = EigenvalueRange::all());

  TridiagonalEigenSolver& computeEigenvectors() {
    eigen_assert(m_isInitialized && "TridiagonalEigenSolver is not initialized.");
    eigen_assert(m_info == Success && "computeEigenvectors() requires a preceding successful eigenvalue computation");
    computeEigenvectorsImpl();
    return *this;
  }

  template <typename DiagType, typename SubdiagType, typename EivalsType>
  TridiagonalEigenSolver& computeEigenvectors(const MatrixBase<DiagType>& diag, const MatrixBase<SubdiagType>& subdiag,
                                              const MatrixBase<EivalsType>& eigenvalues);

  const VectorType& eigenvalues() const {
    eigen_assert(m_isInitialized && "TridiagonalEigenSolver is not initialized.");
    return m_eivalues;
  }

  const MatrixType& eigenvectors() const {
    eigen_assert(m_isInitialized && "TridiagonalEigenSolver is not initialized.");
    eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed.");
    return m_eivec;
  }

  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "TridiagonalEigenSolver is not initialized.");
    return m_info;
  }
 
 protected:
  // Sturm counts, bisection targets, and inverse-iteration thresholds are carried in the scalar type
  // itself, so scalars narrower than float (half, bfloat16) compute in float internally: their
  // epsilon cannot represent consecutive integer counts beyond ~1/eps (a 257-element bfloat16
  // identity would report an eigenvalue of "5.19"), and intermediates such as the coincident-shift
  // perturbation quantize away. Results are rounded back to Scalar on output.
  using ComputeScalar = std::conditional_t<(NumTraits<Scalar>::digits() < NumTraits<float>::digits()), float, Scalar>;
 
  void computeEigenvectorsImpl();
 
  VectorType m_eivalues;
  MatrixType m_eivec;
  // The tridiagonal retained by the last compute step, so the staged computeEigenvectors() can
  // re-form T - lambda*I without the caller having to keep diag/subdiag alive.
  VectorType m_diag;
  VectorType m_subdiag;
  // Narrow scalars only: the unrounded float eigenvalues of the last computeEigenvalues(), kept so
  // the staged eigenvector pass shifts by them rather than by the Scalar-rounded m_eivalues --
  // shifts quantized to Scalar collapse neighbouring eigenvalues and yield near-parallel vectors.
  // Empty whenever m_eivalues did not come from this solver's own bisection.
  Matrix<ComputeScalar, Dynamic, 1> m_eivaluesc;
  ComputationInfo m_info = InvalidInput;
  bool m_isInitialized = false;
  bool m_eigenvectorsOk = false;
};
 
template <typename Scalar_>
template <typename DiagType, typename SubdiagType>
TridiagonalEigenSolver<Scalar_>& TridiagonalEigenSolver<Scalar_>::computeEigenvalues(
    const MatrixBase<DiagType>& diag, const MatrixBase<SubdiagType>& subdiag, const EigenvalueRange& range) {
  static_assert(internal::is_same<typename DiagType::Scalar, Scalar>::value &&
                    internal::is_same<typename SubdiagType::Scalar, Scalar>::value,
                "diag and subdiag must have the solver's scalar type");
  const Index n = diag.size();
  EIGEN_UNUSED_VARIABLE(n);
  eigen_assert(subdiag.size() == (n > 0 ? n - 1 : 0) && "sub-diagonal must have one fewer entry than the diagonal");
 
  m_diag = diag;
  m_subdiag = subdiag;
  m_eigenvectorsOk = false;
 
  // Reject non-finite input up front (cf. SelfAdjointEigenSolver::computeFromTridiagonal()): the
  // Gershgorin bracketing below would otherwise iterate on NaN brackets and return garbage with
  // info() == Success. allFinite() rather than a max-reduction, whose default PropagateFast
  // semantics may not surface a NaN.
  if (!(m_diag.allFinite() && m_subdiag.allFinite())) {
    m_eivalues.resize(0);
    m_eivaluesc.resize(0);
    m_info = NoConvergence;
    m_isInitialized = true;
    return *this;
  }
 
  if (internal::is_same<Scalar, ComputeScalar>::value) {
    internal::tridiagonal_bisection(m_diag, m_subdiag, range, RealScalar(0), m_eivalues);
    m_eivaluesc.resize(0);
  } else {
    // Narrow scalar: bisect in float and round the eigenvalues back (see ComputeScalar). The
    // unrounded values are retained for the staged eigenvector pass (see m_eivaluesc).
    const Matrix<ComputeScalar, Dynamic, 1> cdiag = m_diag.template cast<ComputeScalar>();
    const Matrix<ComputeScalar, Dynamic, 1> csubdiag = m_subdiag.template cast<ComputeScalar>();
    internal::tridiagonal_bisection(cdiag, csubdiag, range, ComputeScalar(0), m_eivaluesc);
    m_eivalues = m_eivaluesc.template cast<Scalar>();
  }
  m_info = Success;
  m_isInitialized = true;
  return *this;
}
 
template <typename Scalar_>
template <typename DiagType, typename SubdiagType, typename EivalsType>
TridiagonalEigenSolver<Scalar_>& TridiagonalEigenSolver<Scalar_>::computeEigenvectors(
    const MatrixBase<DiagType>& diag, const MatrixBase<SubdiagType>& subdiag,
    const MatrixBase<EivalsType>& eigenvalues) {
  static_assert(internal::is_same<typename DiagType::Scalar, Scalar>::value &&
                    internal::is_same<typename SubdiagType::Scalar, Scalar>::value &&
                    internal::is_same<typename EivalsType::Scalar, Scalar>::value,
                "diag, subdiag and eigenvalues must have the solver's scalar type");
  const Index n = diag.size();
  EIGEN_UNUSED_VARIABLE(n);
  eigen_assert(subdiag.size() == (n > 0 ? n - 1 : 0) && "sub-diagonal must have one fewer entry than the diagonal");
  eigen_assert(eigenvalues.size() <= n && "cannot request more eigenvectors than the size of the matrix");
 
  m_diag = diag;
  m_subdiag = subdiag;
  m_eivalues = eigenvalues;
  m_eivaluesc.resize(0);  // the caller's Scalar eigenvalues are authoritative here
 
  // Reject non-finite input up front, mirroring computeEigenvalues(): inverse iteration on NaN/Inf
  // data (or NaN shifts) would return NaN vectors with info() == Success.
  if (!(m_diag.allFinite() && m_subdiag.allFinite() && m_eivalues.allFinite())) {
    m_eivalues.resize(0);
    m_eivec.resize(m_diag.size(), 0);
    m_info = NoConvergence;
    m_isInitialized = true;
    m_eigenvectorsOk = false;
    return *this;
  }
 
  computeEigenvectorsImpl();
  return *this;
}
 
template <typename Scalar_>
void TridiagonalEigenSolver<Scalar_>::computeEigenvectorsImpl() {
  const Index n = m_diag.size();
  const Index m = m_eivalues.size();
  m_eivec.resize(n, m);
  Index nonconv;
  if (internal::is_same<Scalar, ComputeScalar>::value) {
    nonconv = internal::tridiagonal_inverse_iteration(m_diag, m_subdiag, m_eivalues, m_eivec);
    // Refine the eigenvectors of any genuinely degenerate cluster (Rayleigh-Ritz). The eigenvalues
    // are left unchanged, and a non-degenerate spectrum is untouched.
    internal::tridiagonal_rayleigh_ritz_refine(m_diag, m_subdiag, m_eivalues, m_eivec);
  } else {
    // Narrow scalar: iterate and refine in float, then round the vectors back (see ComputeScalar).
    // Shift by the retained unrounded eigenvalues when they exist (the staged path): shifts rounded
    // to Scalar collapse neighbouring eigenvalues and produce near-parallel vectors.
    const Matrix<ComputeScalar, Dynamic, 1> cdiag = m_diag.template cast<ComputeScalar>();
    const Matrix<ComputeScalar, Dynamic, 1> csubdiag = m_subdiag.template cast<ComputeScalar>();
    const Matrix<ComputeScalar, Dynamic, 1> cw =
        (m_eivaluesc.size() == m) ? m_eivaluesc : m_eivalues.template cast<ComputeScalar>();
    Matrix<ComputeScalar, Dynamic, Dynamic> cvec(n, m);
    nonconv = internal::tridiagonal_inverse_iteration(cdiag, csubdiag, cw, cvec);
    internal::tridiagonal_rayleigh_ritz_refine(cdiag, csubdiag, cw, cvec);
    m_eivec = cvec.template cast<Scalar>();
  }
 
  // Like LAPACK xSTEIN, report NoConvergence if any eigenvector failed the inverse-iteration growth
  // test (the columns are still returned, best effort, so eigenvectors() remains usable).
  m_info = (nonconv == 0) ? Success : NoConvergence;
  m_isInitialized = true;
  m_eigenvectorsOk = true;
}
 
}  // namespace Eigen
 
#endif  // EIGEN_TRIDIAGONAL_EIGENSOLVER_H