Non-recursive solution of sparse block Hessenberg systems

A double-phase algorithm, based on the block recursive LU decompo- sition, has been recently proposed to solve block Hessenberg systems with sparsity properties. In the first phase the information related to the fac- torization of A and required to solve the system, is computed and stored. The solution of the system is then computed in the second phase. In the present paper the algorithm is modified: the two phases are merged into a one-phase version having the same computational cost and allowing a saving of storage. Moreover, the corresponding non recursive version of the new algorithm is presented, which is especially suitable to solve infi- nite systems where the coefficient matrix dimension is not a-priori fixed and a subsequent size enlargement technique is used. A special version of the algorithm, well suited to deal with block Hessenberg matrices having also a block band structure, is presented. Theoretical asymptotic bounds on the computational costs are proved. They indicate that, under suit- able sparsity conditions, the proposed algorithms outperform Gaussian elimination. Numerical experiments have been carried out, showing the- effectiveness of the algorithms when the size of the system is of practical interest.