Abstract. We observe that the cyclic reduction algorithm leaves unchanged the structure of a block Toeplitz matrix in block Hessenberg form. Based on this fact, we devise a fast algorithm for computing the probability invariant vector of stochastic matrices of a wide class of Toeplitz-like matrices arising in queueing problems. We prove that for any block Toeplitz matrix $H$ in block Hessenberg form it is possible to carry out the cyclic reduction algorithm with $O(k^3n+k^2n\log n)$ arithmetic operations, where $k$ is the size of the blocks and $n$ is the number of blocks in each row and column of $H$. The probability invariant vector is computed within the same cost. This substantially improves the $O(k^3n^2)$ arithmetic cost of the known methods based on Gaussian elimination. The algorithm, based on the FFT, is numerically weakly stable. In the case of semi-infinite matrices the cyclic reduction algorithm is rephrased in functional form by means of the concept of generating function and a convergence result is proved. An implementation of the algorithm for the infinite case is presented together with the results of the numerical experiments which confirm the effectiveness of our approach.

<postscript>  files of the paper.