title: An accurate and time-parallel rational exponential integrator for hyperbolic and oscillatory PDEs
authors: Marco Caliari Lukas Einkemmer Alexander Moriggl Alexander Ostermann
text: 
  
   
  [t1]The numerical computations are performed using the GPU cluster GPU3 at the University of Innsbruck.
  
  1]Marco Caliarimarco.caliari@univr.it
  
  2]Lukas Einkemmerlukas.einkemmer@uibk.ac.at
  
  2]Alexander Morigglalexander.moriggl@uibk.ac.at
  
  2]Alexander Ostermannalexander.ostermann@uibk.ac.at
  
  
  [1]Department of Computer Science, University of Verona, Italy
  [2]Department of Mathematics, University of Innsbruck, Austria
  
  
  
  
  
  
      Rational exponential integrators (REXI) are a class of numerical methods that are well suited for the time integration of linear partial differential equations with imaginary eigenvalues. Since these methods can be parallelized in time (in addition to the spatial parallelization that is commonly performed) they are well suited to exploit modern high performance computing systems. In this paper, we propose a novel REXI scheme that drastically improves accuracy and efficiency. The chosen approach will also allow us to easily determine how many terms are required in the approximation in order to obtain accurate results. We provide comparative numerical simulations for a shallow water equation that highlight the efficiency of our approach and demonstrate that REXI schemes can be efficiently implemented on graphic processing units.
  
  
  rational exponential integratorsparallel in timehyperbolic problemshighly oscillatory problemsGPU computing
  
  An accurate and time-parallel rational exponential integrator for hyperbolic and oscillatory PDEst1
      [
      December 30, 2023
  ===================================================================================================
  
  
  
  
  
  § INTRODUCTION
  
  
  
  In this work, we are interested in simulating linear partial differential equations (PDEs) with purely imaginary eigenvalues of large modulus (i.e. stiff problems). That is, we consider 
  
      ∂_t f = Af,    f(0) = f_0,   σ(A) ⊂ iℝ,
  
  where σ(A) is the spectrum of A. Such problems arise, for example, in quantum dynamics (e.g. the Schrödinger equation) and wave propagation (e.g. the Helmholtz equation). In addition, solving these problems is an integral part of applying exponential integrators or splitting methods to a large number of nonlinear problems ranging from plasma physics to electrodynamics.
  
  Due to the stiff nature of these equations, explicit time stepping methods are forced to take excessively small time steps in order to remain stable. Thus, implicit schemes (e.g. <cit.>), implicit-explicit IMEX schemes (e.g. <cit.>), exponential integrators (e.g. <cit.>), or splitting methods (e.g. <cit.>) are commonly used. These methods enjoy better stability properties and can thus, in principle, take large time steps. However, for highly oscillatory problems, the maximal time step size of implicit methods is still severely limited by the fact that such methods need to resolve the oscillations.
  
  Approximately a decade ago, the hardware used to run such simulations has undergone a paradigm shift. Due to the fact that frequency scaling has essentially ended at that point, the main way to increase performance has been to add more parallelism. Desktop computers now routinely have 16 cores and large vectorization units. This trend is even more pronounced in high performance computing systems, where supercomputers with millions of threads are now in operation. In addition, graphic processing units (GPUs) have come to the forefront as they are able to outperform central processing units (CPUs) for many scientific computing tasks. This advantage is achieved by providing a massive parallel system. In fact, a sequential program on a GPU would be slower than on a CPU. Therefore, methods which are highly parallelizable and work well on these new computer architectures are needed to take advantage of their computational power. A significant body of research has been accumulated in recent years that considers numerical methods that are well suited for such systems (see, e.g., <cit.>). More specifically, in the context of exponential integrators we refer to <cit.>.
  
  Since time stepping methods are inherently sequential, they generally can not be parallelized in time (although some parallelism can be extracted by constructing methods with parallel stages; see, e.g., <cit.>). Even exponential integrators of high order, which are able to perform large time steps, are mainly implemented sequentially (in time). For example, polynomial approximations, as in <cit.>, require us to calculate a sequence of matrix-vector products which can not be done in parallel. A similar argument holds for Krylov approximations <cit.>. Of course, these methods are parallelizable in space. But in some situations the scalability is limited and, at some point, increasing the number of computing cores does not further reduce the simulation time <cit.>. Therefore, such an approach can not fully exploit modern computer hardware.
  
  A novel idea to overcome this problem are so-called Rational Exponential Integrators (REXI) schemes, which were introduced in <cit.>. The basic idea of these methods is to approximate e^tA by a linear combination of simple rational functions. The advantage of REXI methods is that the corresponding terms can be calculated independently of each other. Therefore, these methods are highly parallelizable in time. It is worth mentioning that REXI is markedly distinct from time parallelization schemes such as the parareal and similar methods (see, e.g., <cit.>). In the latter case, a coarse and a fine time integrator are combined to achieve parallelism within an iterative procedure, while in the former the action of the matrix exponential is directly approximated in a way that is amendable to parallelization.
  
  In this paper our goal is twofold. First, we propose a modification to the original REXI scheme that drastically improves the accuracy and efficiency of the method (section <ref>). The proposed method also allows us to easily determine how many terms the approximation requires in order to obtain accurate results. These theoretical considerations are then confirmed by numerical experiments in section <ref>. Second, we demonstrate that these types of methods can be efficiently implemented on massively parallel computer architectures. Specifically, we demonstrate an implementation on modern GPUs that yields a drastic speedup compared to the corresponding CPU implementation for the shallow water equations (section <ref>).
  
  
  
  
  § THE ORIGINAL AND IMPROVED REXI SCHEMES
  
  
  
  In this section we discuss the derivation of the REXI schemes. Moreover, we reveal some problems of the original scheme in the matrix case and show how they can be eliminated with our new formulation.
  
  
  
   §.§ The scalar case
  
  
  In this section we will give a brief summary of how REXI approximates e^ix with x∈ℝ. For more details, see <cit.>. Our notation is the same as in <cit.>. The three main steps are as follows:
  
  0em 
  
    * Approximate e^ix by a sum of Gaussian functions.
  
    * Approximate each Gaussian function by a sum of rational functions.
  
    * Combine 1. and 2. to approximate e^ix by a sum of rational functions.
  
  
  
  
    §.§.§ Step 1
  
  
  We start by writing e^ix as a linear combination of Gaussian functions <cit.>
  
      e^ix + ϵ = f(x) = ∑_m = -∞^∞ b_m ψ_h(x + mh),
   
  where
  
      ψ_h(x) = 1/√(4π)exp(-x^2/4h^2)
  
  and ϵ = e^h^2∑_k ≠ 0e^-4π^2k^2, see Appendix B. Since ϵ, in general, is smaller than machine precision (see Appendix B or <cit.>), we neglect it in the following.
  
  The parameter h defines the numerical support of the Gaussian function. Our goal now is to determine the coefficients b_m. To that end equation (<ref>) is transformed to Fourier space. Using the shift property of the Fourier transform, it follows that
  
      f̂(ω)/ψ̂_h(ω) = ∑_m = -∞^∞ b_m e^2π i mhω,
  
  and therefore we obtain 
  
      b_m = h∫_-1/2h^1/2h e^-2π i m h ωf̂(ω)/ψ̂_h(ω)  dω.
  
  Since f(x) = e^ix, we have f̂(ω) = δ(ω - 1/2π), where δ(x) is the Dirac distribution. Moreover 
  
      ψ̂_h(ω) = he^-4π^2h^2ω^2
  
  and therefore the sought after coefficients are given by
  
      b_m = e^-imhe^h^2.
  
  Here, the first constraint on h arises. If h > π then the integral in (<ref>) is 0. In practice, this parameter has to be even smaller to produce accurate approximations. Its value is discussed in Section <ref>.
  
  Finally, the sum in (<ref>) has to be truncated:
  
      e^ix≈∑_m = -M^M b_m ψ_h(x + mh).
   
  Therefore, REXI depends on two parameters: M and h. As mentioned earlier, the parameter h < π defines the numerical support of (<ref>) and the parameter M controls, together with h, the interval on which the approximation is sought. It can be shown (<cit.> and Appendix B), that (<ref>) produces an accurate approximation if 
  
      |x| ≤ (M-11)h.
  
  
  
  
    §.§.§ Step 2
  
  
  The second step consists in approximating the Gaussian function ψ_h(x) as a sum of rational functions:
  
  
      ψ_h(x) ≈(∑_l = -L^La_l/ix/h + μ + il) =: R(x/h),
  
  where μ∈ℝ and a_l ∈ℂ are coefficients. To determine these parameters, the authors in <cit.> first approximate ψ̂_1(ω) by a linear combination of exponential functions, 
  
  
      ψ̂_1(ω) ≈∑_j=1^J b_j e^θ_j ω,
  
  using the Adamyan–Arov–Krein theory (see <cit.> for more details). Moving back to physical space, they then obtain
  
      ℱ^-1(∑_j=1^J b_j e^θ_j ω) = -2(∑_j=1^Jb_j/2π ix+θ_j) ≈ψ_1(x),
  
  where (θ_j) < 0. In the original  scheme <cit.> only the coefficients θ_j are computed. Then, they set μ = min_j (θ_j / 2π) and then look for an approximation of ψ_1(x) of the form (<ref>) with h=1,
  
  
      R(x) = (∑_l = -L^L a_l/ix + μ + il),
  
  where L = 11 is chosen. This slight modification of (<ref>) gives the advantage that all the shifted Gaussians share the same poles, which reduces the terms in the next step. Then, they determine the coefficients a_l by minimizing the following l^∞ error
  
      max_x_k|1/√(4π)e^-x_k^2/4 - R(x_k)|
  
  using a set of points x_k ∈ [-30,30]. In Table 1 of <cit.>  the values of the coefficients a_l and μ are given. The corresponding approximation error is less than 7 × 10^-13. 
  
  We now propose a strategy to reduce the approximation error by using a different approach to determine the coefficients. First, we observe that the function ψ_1(x) is a symmetric function, i.e. ψ_1(x)=ψ_1(-x). Therefore, the same property should hold true for the corresponding approximation R(x), which is not satisfied in <cit.>. We remark that R(x) is symmetric if a_l = a_-l for all l. To enforce that this property holds also numerically, we rewrite (<ref>) as
  
      R(x) = a_0μ/x^2 + μ^2 + ∑_l = 1^L 2μ(a_l)
      (μ^2 + l^2 + x^2) + 2l(a_l)(μ^2 + l^2 - x^2)/x^4 + 2(μ^2 - l^2)x^2 + (μ^2 + l^2)^2.
  
  
  Then, to find the coefficients a_l we minimize the following l^2 error
  
      ∑_k=1^K (1/√(4π)e^-x_k^2/4 - R(x_k))^2
  
  on K points x_1, … , x_K. 
  Since R(x) is a linear combination of the coefficients, the approximation can be written as
  
      R(x_k) = G(x_k,μ,L)y ≈ψ_1(x_k),
  
  where G(x_k,μ,L) ∈ℝ^1 × (2L+1) and y ∈ℝ^2L+1 is of the form
  
      y = [a_0, (a_1), (a_2), …, (a_L), (a_1), …, (a_L)]^T.
  
  
  Note that A(x_k,μ,L) is defined by (<ref>).
  Extending this notation, we obtain a linear system 
  
      A(X_K,μ,L)y = ψ_1(X_K)
  
  where 
  
      X_K = [x_1, …, x_K]^T,
  
  
      A(X_K,μ,L) ∈ℝ^K×(2L+1),
  
  
      A(X_K,μ,L) = 
      [ A(x_1,μ,L);          ⋮; A(x_K,μ,L);            ],
  
  
      ψ_1(X_K) = [ψ_1(x_1), …, ψ_1(x_K)]^T ∈ℝ^K.
  
  
  The reason why we chose to minimize (<ref>) instead of (<ref>) is because it can be solved easier and faster. To minimize (<ref>) we compute the least square solution of the corresponding linear system, where the points x_k are calculated iteratively. We start with x_1 = 0 and for selecting the next point x_k+1 we use the same strategy that is used for minimizing the error in interpolation with Leja points <cit.>. For L = 24 we obtained the coefficients listed in Table <ref>. The coefficient μ is determined such that a high accuracy is obtained. For this choice the error in the maximum norm is less than 8× 10^-15. 
  
  
  
  
  
    §.§.§ Step 3
  
  
  The third step is the combination of steps 1 and 2:
  
      e^ix≈∑_m = -M^M b_m ∑_l = -L^L (
      ha_l/ix + h(μ + i(m+l))).
  
  
  For computational efficiency equation (<ref>) should be rewritten as a single sum. This can be done in different ways. A possibility is to split up the coefficients b_m into their real and imaginary parts and pull them inside the approximation of the Gaussian function. This is done in <cit.>. Let  n = m + l, N = M + L, α_n = h(μ + in) and set
  
      β_n^Re = h ∑_k = L_1(n)^L_2(n) a_k(b_n-k)
  
  and
  
      β_n^Im = h ∑_k = L_1(n)^L_2(n) a_k(b_n-k),
   
  where L_1(n) = max(-L,n-M), L_2(n) = min(L,n+M). This leads to the following form 
  
  
      (ix) := ∑_n = -N^N(β_n^Re/ix + α_n) + i(β_n^Im/ix + α_n),
   
  which serves as a definition of the numerical approximation (ix) to e^ix.
  
  Another possibility to rewrite (<ref>) is to compute the real part of the approximation of the Gaussian function. In contrast to the previous reformulation, we make now explicit use of the fact that x has to be real. We set
  
      c_1,n = h ∑_k = L_1(n)^L_2(n)(a_k)b_n-k
   
  and
  
      c_2,n = h ∑_k = L_1(n)^L_2(n)(a_k)b_n-k,
   
  and obtain an equivalent form of (<ref>) which we call , thus
  
      (ix) := ∑_n = -N^Nc_1,nhμ + c_2,n(x+hn)/(α_-n-ix)(α_n +ix).
  
  Recall that N = M+L and α_n = h(μ+in).  In the scalar case both simplifications take roughly the same computational effort and the results are equivalent for real x. This, however, is not true in the matrix case (as we will see in the next section).
  
  We now explain how (<ref>) can be used to compute e^τ A for a given square matrix with purely imaginary eigenvalues. Assume that A is diagonalizable by a matrix V, i.e.,
  
      A = VEV^-1
  
  where E is a diagonal matrix. Then it follows that e^τ A = Ve^τ EV^-1. The diagonal entries of e^τ E can be computed componentwise by REXII
  
      e^τ E = 
        [       ⋱                ;         e^iλ_jτ;                       ⋱ ]≈[       ⋱                ;         (iλ_jτ);                       ⋱ ]
      =: (τ E),
  
  where iλ_j are the eigenvalues of A.  Finally we obtain
  
      e^τ A≈ V(τ E)V^-1.
  
  This is a well known technique to extend scalar functions to matrices, where the scalar functions are applied to the spectrum of the matrix. Therefore, τ A can be substituted for ix in (<ref>) since A is a matrix with purely imaginary eigenvalues. 
  
  
  
   §.§ The new scheme  for matrices
  
  
  
  Let A be a matrix with purely imaginary eigenvalues. As explained before, we can substitute τ A for ix in (<ref>). This yields 
  
      e^τ A≈∑_n = -N^N (c_1,nhμ I + c_2,n(-iτ A+hnI))(α_-nI-τ A)^-1(α_n I+τ A)^-1.
  
  
  
  This scheme can be made more efficient by defining C_1,n = (c_1,nhμ + c_2,nhn) and C_2,n = ic_2,n. As
  
      (C_1,nI - C_2,nτ A)   (α_-nI - τ A)^-1 = 
         C_2,n((C_1,n/C_2,n - α_-n)(α_-nI - τ A)^-1 +I )
   
  we obtain the following extension of  for matrices with purely imaginary eigenvalues
  
      (τ A) := ∑_n = -N^N C_2,n((C_1,n/C_2,n - α_-n)(α_-nI - τ A)^-1 +I )(α_n I+τ A)^-1.
  
  Recall that N = M+L where L=24 and α_n = h(μ + in).  We are now in a position to show how for  the accuracy of the scalar case translates to the matrix case. This is the content of Theorem 1. Its proof is a consequence of (<ref>) and (<ref>).
  
  
  
  
  Let A be a square matrix and suppose that σ (A) ⊂ iℝ. If A is diagonalizable, i.e., A = VEV^-1 with E = (iλ_j), then
  
      (τ A) - e^τ A_∞≤(V)·max_λ_j ∈σ(A)|(iτλ_j) - e^iτλ_j|
  
  where (V) = V_∞·V^-1_∞. 
  
  
  
  Regarding the error in the matrix case, by Theorem <ref> we have to estimate  in the scalar case for every eigenvalue iτλ_j of τ A. In the scalar case  is accurate if condition (<ref>) holds for every iτλ_j, namely |τλ_j| ≤ (M-11)h. This implies in the matrix case, that if we choose M and h such that
  
      τρ(A) ≤ (M-11)h
  
  holds true, where ρ(A) is the spectral radius of A, we are guaranteed to obtain results close to machine precision. In practice, since  in general is applied to a vector, also the error from solving the linear systems has to be taken into account.
  
  
  From (<ref>) we observe that it might be possible to reduce the amount of work by using a shifted matrix, i.e. A' = A-ν I.  If ρ(A') < ρ(A) then less terms are needed to obtain an accurate approximation of e^A'. Note that the matrix e^A can easily be recovered from e^A' as follows
  
      e^A = e^A -ν I + ν I = e^ν e^A - ν I = e^ν e^A'.
  
  Now, since in our case all eigenvalues are purely imaginary, we have iλ_j ∈ i[ζ_1,ζ_2]. Thus to reduce ρ(A) we use the shift ν = iζ_1 + ζ_2/2. This implies that ρ(A') = ζ_1-ν and iλ'_j ∈ i[ζ_1-ν,ζ_2-ν]. For skew symmetric real matrices, it is not necessary to perform a shift since all eigenvalues arise in complex conjugate pairs, and therefore ζ_1 = -ζ_2.
  
  
  
  If the matrix A = VEV^-1 is skew Hermitian, then the matrix V is unitary. If moreover, the error in (<ref>) is estimated in 2-norm, then (V) = 1 and therefore the error in the matrix case is the same as in the scalar case for the eigenvalues of A.
  
  
  
  Since  in general is used to evaluate the action of the matrix exponential applied to a vector f_0, the main cost is to solve two linear systems for each summation term. We are able to reduce the cost if the entries of A and f_0 are real. In particular, we observe that c_1,n = c_1,-n, c_2,n = -c_2,-n and α_n = α_-n. Therefore, we obtain
  
      e^τ A≈( ∑_n = 0^N Γ_nC_2,n((C_1,n/C_2,n - α_-n)(α_-nI - τ A)^-1 +I )(α_n I+τ A)^-1),
  
  where Γ_0 = 1 and Γ_n = 2 for 1 ≤ n ≤ N.
  
  
  
  
   §.§ The original  scheme and the differences to 
  
  
  The original REXI scheme was developed for matrices A with real entries. It is based on (<ref>) where τ A is substituted for ix. A further simplification comes from the fact that e^τ A is real which suggests to neglect the imaginary part of (<ref>). The scheme is thus defined as follows
  
      (τ A) = ∑_n = -N^N (β_n^Re(τ A + α_n I)^-1) ).
  
  
  
  
  At first look this approximation to e^τ A might be the preferred one since if  is applied to a real vector f_0, only one linear system has to be solved for each summation term. In contrast,  has two linear systems to solve for each term. But  has some drawbacks with respect to .
  
  
    * First, if A is real with purely imaginary eigenvalues, then V has to be complex and therefore in (<ref>) the matrix V can not be pulled inside the real part of the approximation. Thus, Theorem <ref> does not hold for  and condition (<ref>) does not apply. To improve the accuracy, M has to be increased. The rate of convergence of  can be slow and an extremely large value for M might be  needed if a stringent tolerance is prescribed, as we will show in the numerical experiments.  Therefore, the main advantage of  compared to  is relation (<ref>) that allows us to choose h and M in an appropriate manner to produce the same high accuracy as in the scalar case. The cost of this is that two linear systems have to be solved, and thus the sequential part of the scheme doubles. Despite the increased computational effort, however, our scheme is still significantly faster since we can choose a much smaller M.
  
  
  
    * Second, reducing the sum from 2N+1 terms to N+1 terms, as is done in (<ref>), is not feasible for the original REXI scheme since the coefficients a_l in <cit.> are not exactly equal to a_-l. Doing this, as in <cit.>, result in a reduction of accuracy in the approximation of the Gaussian functions (from 7 × 10^-13 to 4× 10^-8) and since  is at most as accurate as the approximation of the Gaussian function, the overall accuracy is significantly reduced.
  
  
  
  
  
  
  If A = iB, where B is a real diagonalizable matrix, then the transformation matrix V is real. If moreover f_0 is real, the use of (<ref>) for approximating e^τ A is justified, and we end up with the following scheme:
  
      e^τ Af_0 ≈∑_n=-N^n=N(β_n^(τ A + α_nI)^-1f_0) + i(β_n^(τ A + α_nI)^-1f_0)
  
  We call this scheme  Extended (REXIE). Note that this formulation is more efficient than  as only one linear system has to be solved. Matrices of the above form are, for example, purely imaginary skew Hermitian matrices.
  
  
  
  
  
  
  § NUMERICAL EXAMPLES
  
  
  
  In this section we provide numerical examples that confirm the theoretical considerations laid out in the previous section. We start with the scalar case and then advance to the matrix case. We always use the conjugate symmetric coefficients given in Table <ref> for the implementation of . This is a slight modification of the algorithm in <cit.> but reduces cost and improves accuracy, as explained in the previous section.
  
  
  
   §.§ The scalar case
  
  
  
  Note that in the scalar case, the absolut error is the same as the relative error, since |e^ix| = 1. All the results in this sections are calculated sequentially using GNU Octave.
  
  In order to study the approximation, we fix x and plot the error as a function of h and M. The results are shown in Figure <ref>. We observe that after a certain value of M the error drops immediately to a value close to machine precision. The point at which this happens is well predicted by the bound |x| ≤ (M-11)h. Furthermore, we observe that we obtain the best results for h ∈ [0.3,0.6].  Thus it is not recommended to choose h too large or too small. A small value of h leads to a big value of M,  which increases the computational cost. A large value of h results in reduced accuracy.
  
  
  
  
  We have also investigated the error as a function of x, where we fixed h and let M to be the minimum admissible values given from (<ref>). We observe that for this numerical test the accuracy is always close to machine precision.
  
  
  
   §.§ The matrix case
  
  
  In this section, we analyze the behavior of the different algorithms depending on the parameters h and M for two test matrices A_1 and A_2.
  
  0em 
  
    * The matrix A_1 is the second order finite difference approximation of the advection operator ∂_x with periodic boundary conditions in the spatial domain [a,b] = [0,1]. The discretization step is b-a/n-1 = 1/70. Therefore, A_1 is a skew symmetric matrix with eigenvalues λ_j ∈ i[-70,70]. This matrix has complex eigenvectors, thus we apply  (<ref>). With this matrix  has to solve two linear systems.
  
  
    * The matrix A_2 is the second order finite difference approximation of the free Schrö­din­ger operator i∂_xx in the spatial domain [-1,1] with periodic boundary conditions. The discretization step is 1/35. Thus A_2 is a purely imaginary skew Hermitian matrix with eigenvalues λ_j ∈ i[-4900,0]. Therefore, it is convenient to apply a shift of ν = -2450i. This matrix has real eigenvectors, thus we apply  REXIE (<ref>). With this matrix REXIE has to solve only one linear system per summation term.
   
  As vector f_0 in both cases we used the discretization of (2+cos(2π x))^-1. The results in this section are computed in GNU Octave. We measure the relative error in the l^2 norm
  
  
      err = (A)f_0 - qcr expm(A)f_0_2/qcr expm(A)f_0_2,
  
  where qcrexpm is a Padé approximation of the exponential matrix. The results are shown in Figure <ref>. We clearly see that  is much more accurate for the same M compared to the original  scheme. In addition, the bound (<ref>) predicts the behavior of  very well, which is not the case for the  scheme.
  
  
  
  
  
  
  
  
  
  
  
  
  § LINEAR ROTATING SHALLOW WATER EQUATIONS
  
  
  In this section we apply our proposed REXII scheme and the original REXI scheme to the linear rotating shallow water equations (LRSW) <cit.>. This is the same problem that has been investigated in <cit.>. The LRSW are stated as follows
  
      ∂_t f = Af,    f(0) = f_0,
  
  where the linear operator A is defined as follows
  
      A=
      [    0 -∂_x -∂_y; -∂_x    0    1; -∂_y   -1    0;      ].
  
  The sought after function is f = (η, u ,v)^T, where η is the displacement of the surface height, u the velocity in the x direction and v the velocity in the y direction. The simulation domain is the bi-periodic unit square [0,1]^2, thus periodic in the x direction, f(τ,0,y) = f(τ,1,y), and periodic in the y direction, f(τ,x,0) = f(τ,x,1). The grid resolution is D × D = 128 × 128. To apply REXII (<ref>), we have to solve for each term two linear systems and sum up these two calculated solutions. More specifically we have to compute for 0 ≤ n ≤ N
  
  
      (A + α_n I)g_1,n   = f_0 
      
      (α_-n I - A)g_2,n   = g_1,n
      
      g_3,n   = Γ_nC_2,n(g_1,n + (C_1,n/C_2,n - α_-n)g_2,n) 
      and
      e^Af_0    = ∑_n = 0^N ( g_3,n).
  
  To solve the linear systems, the following strategy is applied (see also <cit.>). Taking the second and the third component of (A + α_n I)g_1,n = f_0 leads to the following equation for the velocities
  
      [ α_n   1;  -1 α_n;     ][ u_1,n; v_1,n;       ]
      = 
      [ u_0; v_0;     ]
      +
      ∇η_1,n.
  
  Inverting the 2×2 linear system yields
  
      [ u_1,n; v_1,n;       ]
      = 
      1/κ_n[ α_n  -1;   1 α_n;     ](
      [ u_0; v_0;     ]
      +
      ∇η_1,n),
  
  where κ_n = 1 + α_n^2. The velocities can then be calculated directly if η_1,n is available. From the first component of (A + α_n I)g_1,n = f_0 we obtain
  
      α_n η_1,n- ∇·[ u_1,n; v_1,n;       ]
      = η_0.
  
  Thus, if we plug (<ref>) into (<ref>) we obtain
  
      Δη_1,n - κ_n η_1,n = r_0,n,
  
  where 
  
      r_0,n   = -κ_n/α_nη_0 + 1/α_nζ_0 - δ_0 
      δ   = u_x + v_y 
      ζ   = v_x - u_y.
  
  Therefore, we first solve the Helmholtz problem (<ref>) for η. Then we obtain the velocities by plugging the gradient of η into (<ref>). To efficiently solve the linear system, all computations will be conducted in Fourier space.
  
  The parameter M of REXI depends on the spectral radius of A. We obtain
  
      ρ(A) = √(2 π^2 D^2+1)≈√(2)π D.
  
  
  We will use different initial conditions to illustrate this example. We compare REXII with REXI and the explicit Runge–Kutta time stepping method of  order 4 (RK4). For the results in the following test scenarios, we use the maximum norm:
  
      err = max_r,s|f_num(τ,x_r,y_s) - f_ref(τ,x_r,y_s)|,
  
  where f_num(τ) is the numerical approximation and f_ref(τ) is a reference solution. For REXI and REXII we perform only one single time step of size τ, and (x_r,y_s) are the grid points of the domain. Since it is not clear for REXI how to choose M and h, we fix h = 0.2 as is done in <cit.> and vary M. Before presenting the numerical results we discuss the parallel implementations for both CPU and GPU based systems.
  
  
  
  
  
   §.§ Implementation
  
  
  
  The calculations in this section are conducted on a single GPU (NVIDIA V100) and separately on a CPU (a dual socket Intel Xeon Gold 5118 server with a total of 32 cores). We also performed the calculations on a NVIDIA TitanV. Since the performance of both NVIDIA cards is quite similar we only report the V100 results here. The code is written in C++ and CUDA 10.0 is used to program the GPU. For the GPU code we use CUFFT <cit.> to perform FFTs and to do the reduction sum at the end we use the CUB library <cit.>. For the CPU code we use the FFTW library <cit.>. The code for the CPU is implemented sequentially over the summation index n, but each iteration is parallelized with OpenMP <cit.>. Initially, we used the BLAS implementation found in Intel MKL <cit.> for the CPU code. However, due to the possibility to aggregate many operations, which is a big advantage for memory bound problems, our OpenMP implementation is actually significantly faster. The GPU code parallelizes over the sum in addition to the spatial parallelization. This is done to exploit the massively parallel architecture of modern GPUs.
  
  
  
  
    §.§.§ CPU Implementation
  
  
  The implementation for the CPU is sequential over the summation range 0:N and is therefore performed exactly how it is described in the previous section. We denote by η̂_1,n the solution of the variable η of the first system and η̂_2,n the solution of the variable η of the second system where n is the index in the sum of REXI. Similar notation holds for the other variables. To perform the calculation of these variables, we use OpenMP with 32 threads.
  
  
  The achieved memory bandwidth given in Table <ref> confirms that the implementation performs as expected.
   
  
  
  
  
    §.§.§ GPU Implementation
  
  
  For the GPU implementation we compute the linear systems in parallel. To facilitate the computation of η̂ we create a matrix of size (N+1) × D^2, where each row corresponds to η̂_n for n = 0 … N. We denote this matrix by E1. In a similar way we create two other matrices for û and v̂, which we call U1 and V1, respectively. We store these matrices in column major format, such that we can apply the fast reduction library CUB to perform the sum over n. Therefore, for REXI we need three such matrices. For REXII we need six of these matrices because we have to solve two different linear systems. We denote these matrices by E1, E2, U1, U2, V1, V2. These matrices are very memory intensive. Each entry requires 16 Bytes of memory, since we are using complex double precision arithmetic. Therefore, for REXII for example, we  need approximately 6 · 16 · D^3 ·√(2)πτ /h bytes of memory to store these six matrices, which corresponds to approximately 90 GBytes for τ = 50, h = 0.5 and D = 128. The GPU we are working with has only 12 GBytes of memory. Therefore, we have to divide the total amount of required memory to store the variables in S parts. Thus we are calculating at each step only a part of E1, namely E1_s = E1(sN̅:(s+1)N̅, :) where s = 0… S-1 and N̅ =  (N+1)/S.
  
  There are six kernel functions required. They carry out the following tasks:
  
  0em 
  
    * computation of E1_s
  
    * computation of U1_s and V1_s
  
    * computation of δ and ζ to solve the second linear system
  
    * computation of E2_s
  
    * computation of the final data
  
    * reduction via CUB library
  
  
  The corresponding algorithm is summarized in Algorithm <ref> and the overall achieved performance is listed in Table  <ref>. We clearly see that the GPU implementation outperforms the CPU implementation by a significant margin.
  
  
  
  
  
  
  
  
  
  
  
   §.§ Wave scenario 1
  
  
  The following initial conditions are used for wave scenario 1:
  
      η(0,x,y)    = sin(4π x)cos(2π y) -15cos(4π x)sin(4 π y) 
      
      u(0,x,y)    = cos(8π x)cos(2π y) 
      
      v(0,x,y)    = cos(4π x)cos(4π y)
  
  This is the same problem as considered in <cit.>. Since we solve the problem in Fourier space, these initial functions are extremely convenient. They are exactly representable in Fourier space with very few terms. This gives us a big advantage when we have to choose the parameter M of REXII. In (<ref>) we can consider D = 6, since all higher modes do not contribute to the solution.
  From (<ref>) we deduce that M is
  
      M = ⌈√(2)·π· 6·τ/h⌉ + 11,
  
  where τ is the final time, since we are performing only one time step. We thus include this problem primarly to provide a comparison to the results obtained in <cit.>. We will conduct an investigation with more realistic initial values in the subsequent sections.
  
  
  
  
  From the results reported in Table <ref>, we recognize that REXII outperforms the original REXI scheme and RK4 in both accuracy and execution time by a large margin. REXI is more comparable to the Runge–Kutta time stepping method of order 4 (RK4) for this initial conditions and RK4 can even outperform the original REXI method in some cases.  
  
  For longer integration times (see the results in Table <ref>), both REXI schemes drastically outperform the explicit RK4 method (as we would expect). In addition, we can see that REXII is much more accurate than the original REXI scheme.
  
  Moreover, we observe that for all numerical methods the GPU implementation significantly outperforms the CPU implementation. For REXI and REXII the speedup ranges from approximately a factor of 7 to a factor of 15.
  
  
   
   
   
   
   
  In the following examples we use M calculated with ρ(A) given by (<ref>) even if it is possible to choose a smaller one as in this example. The reason why we are doing this is that in general we can not expect the initial conditions to be that convenient. 
  
  
  
    §.§.§ Wave scenario 2
  
  
  The same type of initial conditions is used as before with the exception that the frequencies are now much larger. We use:
  
      η(0,x,y)    = sin(32π x)cos(16π y) -15cos(32π x)sin(32 π y) 
      
      u(0,x,y)    = cos(64π x)cos(16π y) 
      
      v(0,x,y)    = cos(32π x)cos(32π y)
  
  
  
  
  Also in this case, see Table <ref>, REXII outperforms the other two methods by a large margin. Here both REXII and REXI work much better than RK4. The reason is that the high frequencies force the explicit time stepping method to take extremely small step sizes. In Table <ref> we show that the onset of convergence strongly depends on the parameter M for this problem. This is expected as  in the matrix case works similarly as in the scalar case. Therefore, the outcome can be compared to the results obtained in Figure <ref>.
   
  
  
  
  
  
    §.§.§ Gaussian scenario
  
  
  The following initial conditions are used for the Gaussian scenario:
  
      η(0,x,y)    = exp(-100((x-0.5)^2 + (y-0.5)^2)) 
      
      u(0,x,y)    = 10^-1sin(64π x)sin(16π y)
      
      v(0,x,y)    = 10^-1sin(32π x)sin(32π y)
  
  
  This initial function η, in contrast to the initial functions in the wave scenarios, is not exactly representable in Fourier space. Therefore, in this case we do not have the advantage of a small spectral radius or that a large part of the frequencies are equal to zero. 
  
  The numerical results for final times τ=1 and τ=50 are shown in Tables <ref> and <ref>, respectively. As before, REXII outperforms the original REXI scheme and RK4 significantly in accuracy. In addition, for both REXII and the original REXI scheme the GPU implementation outperforms the CPU implementation by a factor between 7 and 13.
  
  
  
  
  
  
  
  
  
  § CONCLUSION
  
  
  The original REXI scheme is already a good method to compute the action of the matrix exponential parallel in time. The main downside is that it is not very precise. We proposed a modification of the REXI approach that achieves accuracy close to machine precision at similar or, for some problems, even lower computational cost. The strength of the  methods is the fact that they can be easily parallelized in time (in addition to the commonly used parallelization in space). We have demonstrated this by providing an implementation on massively parallel graphic processing units. The GPU implementation shows a drastic speedup compared to the CPU implementation. 
  
  
  
  § ACKNOWLEDGEMENTS
  
  
  This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 847476. The views and opinions expressed herein do not necessarily reflect those of the
  European Commission.
  
  plain
  
  
  
  
  § APPENDIX A
  
  
  
  
  
  
  
  § APPENDIX B
  
  
  In <cit.> the authors give the idea how to determine an upper error bound, which we will analyze here more in depth. The following expression has to be estimated:
  
  
      |e^ix - ∑_m = -M^M b_m ∑_l = -L^L (ha_l/ix+h(μ + i(m+l)))| .
  
  To do so, the following sum is added and subtracted inside the modulus
  
      ∑_m = -M^M b_m ψ_h(x+mh)
  
  and then the triangle inequality is applied. Thus, we end up with the two terms
  
      |e^ix - ∑_m = -M^M b_m ψ_h(x+mh)|
  
  and
  
      | ∑_m = -M^M b_m ∑_l = -L^L (ha_l/ix+h(μ + i(m+l))) - ∑_m = -M^M b_m ψ_h(x+mh)|.
  
  The term (<ref>) can be estimated using (<ref>):
  
      | ∑_m = -M^M b_m   (∑_l = -L^L (ha_l/ix+h(μ + i(m+l))) - ψ_h(x+mh) )| 
         ≤∑_m = -M^M e^h^2|  R(x/h+m) - ψ_h(x+mh) | 
         ≤ e^h^2(2M+1)δ_2,
  
  where δ_2 is the approximation error of the Gaussian function. In our case, δ_2 = 8× 10^-15, see step 2 in Section <ref>. 
  This bound for (<ref>) is not sharp. It depends linearly on the parameter M, but we observed that after doing several numerical experiments this amount does not increase much when M increases.
  
  
  To deduce an error bound for (<ref>), Poisson's summation formula is applied: 
  
      ∑_m = -∞^∞ϕ_h(x + mh) = 
       1/h∑_k = -∞^∞ e^2π i(k/h)xϕ̂_h (k/h).
  
  For the choice  
  
      ϕ_h(x+mh) = e^-i(x+mh)ψ_h(x+mh)
  
  we get
  
      e^-ix∑_m = -∞^∞ e^-imhψ_h(x+mh) 
       = 1/h∑_k = -∞^∞ e^2π i(k/h)xψ̂_h (k/h +1/2π).
  
  This equation is then multiplied by e^ix and the sum on the right-hand side of the equation is split up, leading to
  
  
      ∑_m = -∞^∞ e^-imhψ_h(x+mh) -e^ix1/hψ̂_h(1/2π) 
          = 1/h∑_k ≠ 0 e^2π ix(k/h + 1/2π )ψ̂_h (k/h +1/2π).
  
  This implies 
  
  
  
      |∑_m = -M^M    e^-imhψ_h(x+mh) -e^ixψ̂_h(1/2π)/h| = 
         = | 1/h∑_k ≠ 0 e^2π ix(k/h + 1/2π )ψ̂_h (k/h +1/2π) 
      -∑_|m|>M e^-imhψ_h(x+mh) | 
         ≤1/h∑_k ≠ 0ψ̂_h (k/h +1/2π) + ∑_|m|>Mψ_h(x+mh).
  
  The coefficients b_m = he^-imhψ̂_h(1/2π)^-1 are given by (<ref>) and ψ̂_h(ω) by (<ref>). Thus, if  the inequality is divided by the positive number ψ̂_h(1/2π)/h it follows
  
  
      |∑_m = -M^M    b_m ψ_h(x+mh) -e^ix| ≤
         ≤1/ψ̂_h(1/2π)(∑_k ≠ 0ψ̂_h (k/h +1/2π) + h ∑_|m|>Mψ_h(x+mh)) 
         ≤ e^h^2(1/h∑_k ≠ 0ψ̂_h (k/h) +  ∑_|m|>Mψ_h(x+mh)) 
         = e^h^2(∑_k ≠ 0 e^-4π^2k^2+∑_|m|>Mψ_h(x+mh)):= e^h^2δ_1.
  
  So the final result is the following:
  
      |e^ix - (ix,M,h)| ≤ e^h^2(δ_1 + (2M+1)δ_2).
  
  Let us now analyse δ_1 = ∑_k ≠ 0 e^-4π^2k^2+∑_|m|>Mψ_h(x+mh). The first sum is neglible. Since e^-4π^2k^2≤ e^-4π^2k it follows that 
  
      ∑_k = 1^∞ e^-4π^2k^2 < ∑_k = 1^∞ e^-4π^2k = 1/1-e^-4π^2-1 ≈ 7.15 × 10^-18.
  
  
  The second sum is more interesting regarding the error. It can be shown that for |x| ≤ (M-m_0)h this sum is also neglible, where m_0 is a constant related to how small the first two terms of the sum ∑_|m|>Mψ_h(x+mh) are. An easy calculation shows that 
  
  
      ψ_h(z) ≤tol  if   |z| ≥ 2h√(-log(√(4π) tol)) = ch.
  
  With tol = 10^-16 we have c ≈ 12. Now it is of interest that ψ_h(x + (M+1)h) < tol and ψ_h(x -(M+1)h) < tol such that these two shifted Gaussian functions have no impact on the approximation of e^ix. This implies, by the previous calculation, that (x - (M+1)h) ≤ -ch and (x +(M+1)h) ≥ ch. Thus, |x| ≤ (M+1-c)h and we may define m_0 := c-1 ≈ 11. 
  
  