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
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?
Dense matrix solving using
dgetrs is way faster in LAPACK, as expected.