  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 transformations efficiently even if no fast transform is
available. We illustrate the performance of the new integrator by solving
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.
