  Exponential integrators have been introduced as an efficient alternative to
explicit and implicit methods for integrating large stiff systems of
differential equations. Over the past decades these methods have been studied
theoretically and their performance was evaluated using a range of test
problems. While the results of these investigations showed that exponential
integrators can provide significant computational savings, the research on
validating this hypothesis for large scale systems and understanding what
classes of problems can particularly benefit from the use of the new techniques
is in its initial stages. Resistive magnetohydrodynamic (MHD) modeling is
widely used in studying large scale behavior of laboratory and astrophysical
plasmas. In many problems numerical solution of MHD equations is a challenging
task due to the temporal stiffness of this system in the parameter regimes of
interest. In this paper we evaluate the performance of exponential integrators
on large MHD problems and compare them to a state-of-the-art implicit time
integrator. Both the variable and constant time step exponential methods of
EpiRK-type are used to simulate magnetic reconnection and the Kelvin--Helmholtz
instability in plasma. Performance of these methods, which are part of the EPIC
software package, is compared to the variable time step variable order BDF
scheme included in the CVODE (part of SUNDIALS) library. We study performance
of the methods on parallel architectures and with respect to magnitudes of
important parameters such as Reynolds, Lundquist, and Prandtl numbers. We find
that the exponential integrators provide superior or equal performance in most
circumstances and conclude that further development of exponential methods for
MHD problems is warranted and can lead to significant computational advantages
for large scale stiff systems of differential equations such as MHD.
