  In the present paper we consider numerical methods to solve the discrete
Schr\"odinger equation with a time dependent Hamiltonian (motivated by problems
encountered in the study of spin systems). We will consider both short-range
interactions, which lead to evolution equations involving sparse matrices, and
long-range interactions, which lead to dense matrices. Both of these settings
show very different computational characteristics. We use Magnus integrators
for time integration and employ a framework based on Leja interpolation to
compute the resulting action of the matrix exponential. We consider both
traditional Magnus integrators (which are extensively used for these types of
problems in the literature) as well as the recently developed commutator-free
Magnus integrators and implement them on modern CPU and GPU (graphics
processing unit) based systems.
  We find that GPUs can yield a significant speed-up (up to a factor of $10$ in
the dense case) for these types of problems. In the sparse case GPUs are only
advantageous for large problem sizes and the achieved speed-ups are more
modest. In most cases the commutator-free variant is superior but especially on
the GPU this advantage is rather small. In fact, none of the advantage of
commutator-free methods on GPUs (and on multi-core CPUs) is due to the
elimination of commutators. This has important consequences for the design of
more efficient numerical methods.
