Many important and complex real-world physical phenomena can be effectively described by stiff systems of Advection--Diffusion--Reaction (ADR) equations. Consequently, developing efficient numerical integrators is crucial for enabling fast and reliable simulations. In this context, exponential time integration schemes have attracted significant attention in recent years due to their strong performance in stiff regimes. This talk focuses on a specific class of such methods, i.e., the so-called directionally split exponential integrators, which are well-suited for systems exhibiting a d-dimensional Kronecker sum structure. In particular, we will present efficient strategies for computing the arising actions of exponential-like matrix functions by means of tailored tensor-matrix operations, such as the mu-mode product and the Tucker operator. These techniques are highly scalable on modern computing architectures, including Graphics Processing Units (GPUs). Time permitting, we will briefly touch on the computational aspects related to the GPU implementation. Finally, we demonstrate the effectiveness of the proposed approach through numerical experiments involving the time integration of semidiscretized ADR systems that give rise to Turing patterns.