I am interested in computing the solution of a lage system of ODEs using a krylov method as in . Such method involve functions related to the exponential (the so-called $\varphi$-functions). It essentially consists of computing the action of the matrix function by constructing a Krylov subspace using Arnoldi iteration and projecting the function on this subspace. This reduce the problem to compute the exponential of a much smaller Hessenberg matrix.
I am aware that there are several algorithms to compute the exponential (see  and references therein). I wonder if there is a special algorithm to calculate the exponential that can takes advantage of the fact that the matrix is Hessenberg ?
 Sidje, R. B. (1998). Expokit: a software package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-156.
 Moler, C., & Van Loan, C. (1978). Nineteen dubious ways to compute the exponential of a matrix. SIAM review, 20(4), 801-836.
 Moler, C., & Van Loan, C. (2003). Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM review, 45(1), 3-49.