  In this paper, we propose a $\mu$-mode integrator for computing the solution
of stiff evolution equations. The integrator is based on a $d$-dimensional
splitting approach and uses exact (usually precomputed) one-dimensional matrix
exponentials. We show that the action of the exponentials, i.e. the
corresponding batched matrix-vector products, can be implemented efficiently on
modern computer systems. We further explain how $\mu$-mode products can be used
to compute spectral transforms efficiently even if no fast transform is
available. We illustrate the performance of the new integrator by solving,
among the others, three-dimensional linear and nonlinear Schr\"odinger
equations, and we show that the $\mu$-mode integrator can significantly
outperform numerical methods well established in the field. We also discuss how
to efficiently implement this integrator on both multi-core CPUs and GPUs.
Finally, the numerical experiments show that using GPUs results in performance
improvements between a factor of $10$ and $20$, depending on the problem.
