Skip to content

num::hessenberg_shifted_solve

Defined in header <linear/factorization/hessenberg.hpp>

template <std::floating_point T, class Index>
void hessenberg_shifted_solve(std::complex<T> *y, const T *H, std::complex<T> shift,
                              const std::complex<T> *b, idx n, std::complex<T> *work,
                              Index *piv) noexcept;  // (1)

void hessenberg_shifted_solve(const mat<real> &H, cplx shift, const vec<cplx> &b_tilde,
                              vec<cplx> &y, array<cplx> &M_buf, array<idx> &pivots);  // (2)

(1) Solve \((sI - H)\,y = b\) for a single right-hand side.

(2) Solve the shifted Hessenberg system \((sI - H)y = \tilde b\) in \(O(n^2)\).

Gaussian elimination with partial pivoting on an upper Hessenberg matrix touches one subdiagonal per column, so the factorization is \(O(n^2)\) rather than \(O(n^3)\). This is what makes a resolvent \((sI-A)^{-1}b\) cheap once \(A\) has been reduced once: every subsequent shift reuses the same \(H\).

Parameters

  • H: Upper Hessenberg matrix, n*n.
  • shift: The scalar s.
  • b_tilde: Right-hand side, length n.
  • y: Output, resized to n.
  • M_buf: Scratch, grown to n*n and reusable across calls.
  • pivots: Scratch, grown to n and reusable across calls.