Efficient and stable solution of M-matrix linear systems of (block) Hessenberg form

The problem of solving large linear systems whose coefficient matrix is a sparse M-matrix in block Hessenberg form has recently received much attention, especially for applications in Markov chains and queueing theory. Stewart proposed a recursive algorithm which is shown to be backward stable. Although the theoretical derivation of such an algorithm is very simple, its efficient implementation is logically rather involved. An analysis of its computational cost in the case where the initial coefficient matrix satisfies quite general sparsity properties can be found in [P. Favati et al., Acta Tecnica Acad. Sci. Hungar., 108 (1997-1999), pp. 89-105]. In this paper we devise a different divide-and-conquer strategy for the solution of block Hessenberg linear systems. Our approach follows from a recursive application of the block Gaussian elimination algorithm. For dense matrices, the present method has a computational cost comparable to that of Stewart's algorithm; for sparse matrices it is more efficient and more robust.