  Markov-modulated Brownian motion is a popular tool to model continuous-time
phenomena in a stochastic context. The main quantity of interest is the
invariant density, which satisfies a differential equation associated with the
quadratic matrix polynomial $P(z) = Vz^2-Dz +Q$, where the matrices $V$ and $D$
are diagonal and $Q$ is the transition matrix of a discrete-time Markov chain.
Its solution is typically constructed by computing an invariant pair of $P(z)$
associated with its eigenvalues in the left half-plane, or by solving the
matrix equation $X^2V-XD+Q=0$. We show that these tasks can be solved using a
componentwise accurate algorithm based on Cyclic Reduction, generalizing the
recently appeared algorithms for the linear case ($V=0$). We give a proof of
the numerical stability of our algorithm in the componentwise sense; the same
proof applies to Cyclic Reduction in a more general M-matrix setting which
appears in other applications such as the modelling of QBD processes.
