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

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.

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.