A single server queue with general service time distribution is considered when the input is a Markov modulated Poisson process (MMPP). An algorithmic solution to the transform of the stationary delay and queue length distributions is summarized, and recursive closed-form expressions are obtained for the moments of these distributions. The numerical implementation of these results is discussed in detail with particular reference to an algorithm due to Lucantoni and Ramaswami [11] and its accelerated version due to Ramaswami [19]. This algorithm is shown to be an efficient tool in the matrix-analytic solution of many stochastic models, as various steps for saving considerable amounts of unnecessary computations are identified. A special case of the model where the service time distribution is of phase type is discussed and the stationary queue length distribution at arbitrary times is obtained in matriz-geometric form. Finally, the matrix-geometric and the M/G/I approaches are compared through this special case.