Solving linear system $Ax=b$ with Hessenberg matrix using lapack

Question

I need to solve a linear system of the form $$Ax = b$$ where $A$ is upper Hessenberg matrix with the lower bandwidth equal to 1, $b$ is the RHS vector and $x$ is the solution vector. I have a C++ routine based on "C.B. Moler, Algorithm 423, Linear Equation Solver, C.A.C.M. 15 (1972), p. 274." that solves the system. However, I wish to use LAPACK if available. Since I could not find any LAPACK routines for 'Hessenberg matrix linear system solving', I used the dgetrf and dgetrs routines, but they are way slower than the C++ implementation for Hessenberg matrices. Sample Output is:

Bench-marking Hessenberg matrix solve routines (size = 400)
            Function   Time-taken(s)         Residual
            --------   -------------         --------
          C++ dech()        0.000367               xx
          C++ solh()        0.000140       0.00000004
        LPK dgetrf()        0.004051               xx
        LPK dgetrs()        0.000110       0.00000005

Final results
  C++ factorization is 11.027X faster than LAPACK.
  C++ bf solve is 0.790X slower than LAPACK.

My question is: is there any routine in LAPACK that specializes in handling linear system solves with Hessenberg matrices?

NOTE: Dense matrix solving using dgetrf and dgetrs is way faster in LAPACK, as expected.

Answer 1

There is no (readily available) solver for systems in Hessenberg form inside LAPACK. However, you can build one yourself using LAPACK routines.

Specifically, use Givens rotations to reduce your system to upper triangular form. Then do backward substitution on the transformed system. You construct Given rotations with xROTG. You apply them with xROTF. You solve a triangular system with xTRTRS

There is an inverse iteration for Hessenberg systems built into LAPACK, which could presumably be tricked into solving your system, but I would not bother. Orthogonal transformations will never get you into trouble.

Answer 2

There is an algorithm by Henry ("The Shifted Hessenberg System Solve Computation", 1995) that allows you to combine the Givens rotations and backsubstitution into a single pass, without modifying the matrix (in-place with O(n) storage). We implemented this in the Julia standard LinearAlgebra library and found it to be several times faster than separate Givens QR + xTRTRS solves.