Inexact matrix exponential preconditioner for implicitly restarted Arnoldi method in fluid dynamics stability problems for parallel heterogeneous architecture

We wish to solve the problem of linear stability analysis of stationary solutions for a fluid dynamics problem governed by 3D Navier Stokes equations. The discrete problem is of order O(106) - O(107) so we use parallel computations to perform linear analysis. We aim at the heterogeneous computational architecture based on central (CPU) and graphical (GPU) processing units. In the paper we use the implicitly restarted Arnoldi method (IRAM) to recover the most dangerous eigenvalues. The IRAM has poor convergence without specific transformation of the matrix spectrum. Standard methods require either long computations or inverse of the transformed matrix. This can be troublesome in fluid dynamics problems where matrix is not available explicitly and is non-symmetric, ill-conditioned and can be fully filled. We propose the usage of the inexact matrix exponential with shift-inverse transformation using only matrix-vector product. The inverse is performed using Krylov methods using only matrix-vector product. We prove convergence bounds for the GMRES method and show that this type of matrix transformation accelerates GMRES convergence. We demonstrate application of the method in recovering eigenvalues close to the imaginary axis in some model problems on some matrices from the Matrix Market. We also demonstrate the application of the method for some fluid dynamics problems. We demonstrate convergence, acceleration and efficiency on the multiple GPU cluster.


Introduction
During the bifurcation analysis of fluid dynamics problem we need to solve the eigenvalue problem for the most dangerous part of the spectrum at a certain point in parameter space. We assume that by some procedure we derive the finite dimensional system of equations that is solved on a computer: where v(t) ∈ R n is a vector of a discrete solution of the original system, n is large (typically from 1 · 10 6 to 1 · 10 8 ), M is a linear operator (also called mass matrix, can be an identity E) and F is a discrete operator, that is usually combined of linear and nonlinear parts (such operator F is called Ozeen operators for incompressible flow equations), f ∈ R n is a discrete source term and R is the parameter vector. Boundary conditions are included in M and F operators. For 2 some solution v 0 at R 0 of (1) we derive the linearized system for ininitialismal perturbations y = v − v 0 : where A| v 0 is the linear operator of F at solution v 0 and δR is the variation of parameter vector. Operator A can be represented in operator form (as operator-vector application) or in matrix form, in this case A ∈ R n×n . We would like to know the most dangerous part of the spectrum sp(A) that is responsible for the linear stability loss for δR variation. So have two problems of linear stability. The first one is related to periodic and quasi-periodic solutions and in this case the left hand term is not zero. This problem usually poses no computational challenge for eigenvalue solvers. In case of the zero left hand side we interested in the stability of the stationary solution by finding the part of the spectrum closest to the imaginary axis. This problem poses substational computational challenges for the eigenvalue solvers especially in HPC implementation. In order to solve (1) and (2) with high precision i.e. perform Direct Numerical Simulation (DNS), one must use High Performance Computations (HPC). In this papaer we are aiming at the heterogeneous computational architecture that consits of Central Processing Units (CPUs) and Graphical Processing Units (GPUs) in configuration where number of GPUs and CPUs are installed in a single node and these nodes are connected by fast interconnect. This type of HPC cluster configuration is a standard configuration for almost all Top500 machines. We are mostly interested in GPU dominated machines (up to 8 GPUs installed on one node) because GPUs decrease power consumption and machine cost. Let us give a brief survey of methods that are already developed and successfully used during the linear stability analysis of systems (1), where F is the Ozeen discrete operator. The solution of an eigenvalue problem for (1) in case of incompressible flow is conducted using Krylov subspace methods with exponential estimation preconditioner [20] or shift and inverse preconditioner [1]. One of the first such papers was [7] (and some other papers by these authors), where a 3D Floquet stability analysis was applied to a circular cylinder 2D wake stability problem. The fluid flow solver was a spectral element method for 2D equations where Floquet analysis was used through the introduction of infinitesimal three-dimensional perturbations and derivation of a state-transition periodic monodromy matrix. Such approach can be successfully used for the analysis of first bifurcations but the condition number becomes large still for higher Reynolds numbers. Using this approach bifurcation analysis for incompressible flow problems were successfully conducted in papers, for example in [2,3]. Another successful approach is developed by Sanchez and co-authors, some results can be seen in [4,5,6] also using Krylov methods, namely Arnoldi method. The linear stability of a 3D compressible subsonic flow in a cavity is analyzed in [8] using ARPACK implementation of the implicitly restarted Arnoldi method and a sixth-order compact finite-difference scheme for spatial discretization. Other methods exist that can be successfully applied to these problems. For example the novel method based on two sided Newton method [22,23] is successfully appled to find linear stability and optimal disturbances in many fluid dynamics problems. Such approach can be used for HPC application for systems where a basic Newton-Raphson method is tuned and efficiently implemented.
In all cases we can see that there are no methods of eigenvalue analysis that are specifically aimed on multiple GPU architecture. Our goal is to present a method that is used in many our papers to perform eigenvalue analysis on multiple GPU orientated architecture. For the solution of the eigenvalue problem for (2) we use the use Implicitly Restarted Arnoldi Method (IRAM) [20,21,36]. We chose this method for it's efficiency and need only to have a matrix-vector linear operator to execute Arnoldi process. We solve the following eigenvalue problem: where f () is a function that transforms matrix spectrum in some manner in order to accelerate convergence of the IRAM. If IRAM converges, we obtain diagonal matrix f (Λ) and eigenvector matrix V of the desired part of the spectrum (defined separately for particular f ()). Original eigenvalues are recovered, say, by the Galerkin projection: where desired eigenvalues of the original problem are stored in matrix R. The IRAM is based on the construction of the Krylov subspace (K) dimension p: where u 0 is a starting vector with u 0 = 1 and we designate u n := f (A)u n−1 . Every basis vector in (5) is orthogonalized with respect to all other basis vectors by the Arnoldi process. More information on the IRAM is given in [20] and currently used multiple GPU implementation is described in [36]. The method convergence well if desired part of the spectrum is well separated from the rest of the spectrum and is located on the boundary of the convex hall constructed over the spectrum. In other words, by applying f () we wish to transform the desired part of the spectrum to the boundary and separate it. One usually applies the following transformations: Exponential transformation: Polynomial transformation: Shift-inverse transformation: Cayley transformation: Inverse Kronecker sum transformation: Let us assume that we have two leading real eigenvalues λ 1 and λ 2 close to zero with 0 > λ 1 > λ 2 , |λ 2 | < 1. The matrix exponential transformation (6) would be ideal for the IRAM preconditioning only if it was explicitly available. One would seek transformed eigenvalues with maximum magnitude and recover original eigenvalues with maximum real part. However one usually applies time stepping code to calculate (6) where parameter µ = ∆t is a time step in the numerical method. In this case a error on each iteration is multiplied by the factor if e (λ 2 −λ 1 )∆t ∼ 1, where ∆t << 1. Still this method is used to compute eigenvalues in hydrodynamic stability problems, e.g. [2,3], when unstable eigenvalues are sought. A polynomial transformation (7) applies coefficients c j such that a desired part of the spectrum is increased and the rest is damped. This method is used in ARPACK optionally. However it is difficult to apply this transfirmation to separate eigenvalues close to the origin, so, the method is not used in stability problems that we are interested at.
Shift and inverse (8) and Cayley (9) transformations are, probably, the best candidate to precondition eigenvalues located near a specific point (σ = 0 for those near the origin) for the IRAM. The error will be rediced by the small factor λ 1 /λ 2 for |λ 1 | << 1. One, then seeks eigenvalues with maximum magnitude for transformed spectrum and reovers eigenvalues, closest to point σ. Direct inverse in (8) is never used in practice and one has to solve a system of linear equations in the form (A − σE) u n+1 = u n for the IRAM in order to find a new Krylov vector u n+1 in (5). However, it is well known [9,10] that the condition of A can be very poor, especially for high Reynolds numbers and compressible flows. So the inversion in (8) becomes a huge problem that is to be solved in order to solve the eigenvalue problem. The same is true for (10). The condition of such matrix is worse than that for the Jacobi matrix and problem with preconditioners emerges. This problem is worsten for HPC GPU oriented calculations because GPUs have limited amount of device memory and efficient preconditioners for advection-dominated flows (like ILU factorization) require explicit matrix storage. We present the method that allows us to solve the problem of eigenvalue analysis efficiently on HPC multi GPU computational architecture without explicit matrix storage or matrix graph reconstructions.
The papaer is laid out as follows -first we give the preconditioner theory and theoretical convergence estimates. We also include information on the execution algorithm. Then we describe shift-inverse implementation with specific preconditioner for GPU-oriented computional architecture that we use in our calculations. Next section includes testing problems (matrix market examples and fluid dynamics benchmarks) and computational efficiency estimates of the suggested algorithms. Finaly we present some results for 3D fluid dynamics stability problems that were successfully solved by the suggested algorithm. Next, the conclusion follows.
All calculations are performed using an 8-core (16 in hyper-threading) Intel Xeon E5-2640V2 Ivy Bridge and up to six k40 NVIDIA GPUs assembled in one chassis.

Inexact matrix exponent preconditioning for Arnoldi eigensolver
We denote the left half-plane as C −,0 : {∀(a, b) : a + ib, a ≤ 0} and right half-plane as C + : {∀(a, b) : a + ib, a > 0}. Assume that the spectrum sp(A) ⊂ C −,0 . Then we apply the following two stage transformation for the Krylov basis vector (u n+1 ): The first transformation f 1 (A) = (E + µA) m is the inexact exponent transformation (6) which is executed using explicit Euler's method (or another timestepping method used during the solution of (1)) with m timesteps and µ such that these iterations are stable. The contraction is controled by number of iterations m. Eigenvalues λ 1 > λ 2 are transformed to region (1 − , 0) on the complex plane, where = (1 + ∆tλ 1 ) m . If λ 1 = 0, then it is transformed exactly to (1, 0). By this transformation we have sp(f 1 (A)) ⊂ B (1) (0, 0) (unit 2D ball with the center at the origin). The next transformation is a shift in the manner of shift and inverse transformation. The shift σ defines the part of the spectrum to be sought: σ = 1 − + δ. We are interested in the spectrum near the imaginary axis, so δ is selected such that the solution of the linear system (11) convergence rapidly and the spectrum of the resulting mapping is well separated.
To justify our approach may use the following . Let the GMRES method appled to the solution of Bx = b generate residual r q on the q-th step. Then the GMRES method with q steps convergence to the, at least, given precision Proof. Let us consider the q-th residual r q = b − Bx q of the GMRES method [11]: where P q denotes the set of polynomials of degree at most q and with value one at the origin. Following [11] we derive: International Now, by our assumption sp(A) ⊂ C −,0 so mapping (11) will transfer the whole spectrum of the problem to the B 1 (−1 − δ, 0). Analytical complex function reaches it's maximum at the boundary, i.e. max k |p(λ k )| = max Then we can use trial polynomialĈ q and Chebyshev polynomials C q for the boundary S 1 (−1 − δ, 0) [12] and arrive at: Using (14) and (15) we have: followed by (12) after reduction by 2. This completes the proof.
Technically it means that the convergence doesn't explicitly depend on the size of the matrix A. The dependence on the k(V), as it is said in [13] makes the estimate (13) bad. However the bound of k(V) = M is in the logarithm and logarithmically week and dependence on δ is in the base of the logarithm and very intensive. We demonstrate this fact in figure 1.
Since calculation of matrix exponent is a real operation this method can only be applied to pure real eigenvalue detection. To correct this deficiency we introduce the rotation parameter θ = σ[cos(ϕ) + i sin(ϕ)]. Then the transformation (11) is replaced with the following system: where we have additional basis vectors responsible for real u Re and imaginary u Im mapping.
The angle ϕ ∈ [0, ϕ max ], and maximum angle of rotation ϕ max is known due to the following Proposition 2. Let a timestepper that is used to integrate (1) has no dominant stability limit due to diffusion. Let C be a Courant number from the timestepper and let us use the transformation (11). Then: in transformation (17).
Proof. Let us note that µ is selected such, that µC ≤ 1. If the numerical method is stable, then the spectrum of the operator is bounded by the stability region of the method. Maximum is reached on the imaginary part of eigenvalues since diffusion is not a limiting factor. Knowing the C constant we know a bounding region on an imaginary axis. This bounding region is roughly approximated by a bounding box. The most unstable spectrum is located at the neighbourhood of the imaginary axis. Taking all this into account we consider the boundary eigenvalue λ 1 such, that |Re(λ 1 )| << 1, Re(λ 1 ) < 0, |Re(λ 1 )| << |Im(λ 1 )|. Let us apply the transformation (11) . In this case the location of a transformed eigenvalue f (λ 1 ) on the imaginary axis can be bounded as: So the maximum angle of rotation ϕ max that brings this corner eigenvalue to the real axis by transformation (17) is bounded by (19) from bellow and by π from above.  Figure 1. Illustration of the q vs δ dependence by (12) for a matrix with M = 1 · 10 12 and ε = 1 · 10 11 .
The system increases its size twofold resulting in additional storage space which is a significant drawback. Another drawback results in the introduction of rotational grid that must be used in order to find all leading eigenvalues with nonzero imaginary parts. However this approach, in authors experience, is more stable in practice than the Cayley transformation. The convergence properties don't change since the mapping (17) has no affect on k(V).
In order to recover the original eigenvalues of the problem from the mapped operator (17) we perform the following projections: Desired eigenvalues of the original problem are contained in the matrix R. Notice, that pure real eigenvalues are doubled by this projection.
Let now there exist the part of the spectrum sp 1 (A) : sp(A) ⊃ sp 1 (A) ⊆ C + with leading eigenvalues Re(λ 1 ) > Re(λ 2 ) > 0. It means that the main solution becomes unstable. Then we apply inexact exponential transformation (6) and eigenvalues are transformed as (1 + ∆tλ j ) m , j = {1, 2}. So we controll the separation of the desired spectrum part sp 1 (A) by increasing ∆t or increasing m. This approach for multiple GPU architecture and stability analysis is illustrateded in [36]. For the linear system solution we use restarted GMRES(n) or BiCGstab(L) methods (the latter has better convergence in practice). Schematically, the described algorithm for inexact matrix exponent preconditioner is presented in figure 2.
From this figure we can estimate number of matrix-vector multiplications as: where K IRAM is the number of vectors in Krylov subspace of the IRAM, I Arnoldi is the number of Arnoldi iterations, K GM RES is the number of basis vectors in GMRES, I GM RES is the number of restarts for GMRES method and m is number of iterations in the inexact exponent mapping (11). Typical value of op(M V ) is about 60000-80000 from our experience. This is compareable with the typical M V count for the linear system solution during the standard shift and inverse process with good preconditioning. However, this approach is way more efficient on the multiple  Figure 2. Schematic discription of the inexact exponent shift and inverse IRAM algorithm GPU systems since minimal PCIExpress transfer is used during iterations and only matrix vector applications are used by calling linearised timestepper for the given vector.

Shift-inverse preconditioning on multiple GPUs
We also implemented shift-inverse and Cayley transformations (8) for multiple GPU architecture using simplified approximate preconditioner for the solution of the linear system.
Survey and analysis of preconditioners for fluid dynamics problem or GPU architecture is given in [14,15,16] and for shift inverse methods -in [24,25]. One can observe that the iterative convergence properties of the matrix f (A) are worsten by the decrease of diagonal elements in (8). The preconditioned shift-inverse method is solved iteratively using Krylov-type method (GMRES or BiCGStab(L)). Let x j be an approximation of the eigenvector and B := (A − σE). An iteration of the Krylov method can be formulated with the right preconditioner P as: Olshanskii2012 We construct the preconditioner by using the linearisation on more simple operator with close spectral properties (pseudo-compressibility method for incompressible flows or first order method for compressible flows) with simple Lax-Friedrichs fluxes and simple diffusion discrete operator (7 point Laplace operator in 3D). It is easy to show that the these methods generate diagonally dominant matrices, even for large Reynolds numbers. The Jacobi matrix is found using analytical differentiation of the Lax-Friedrichs flux. Then the matrix can be inverseed by the application of the Block Jacobi or Block Symmetric Gauss-Seidel methods and applied as the preconditioner that gives a rough estimate of the linear operator spectrum. We use Geometric Muligrid (GMG) method [17,18,19,32] to accelerate convergence of the preconditioner. It was tested in many applications on GPU architecture and shown good parallel efficiency.
Let us denote GM G (A, b) as a GMG set of V or W cycles for the solution of the linear system with matrix A and the RHS b. Now we construct the preconditioner. First we split our operator F(u) = N(u, u) + Cu to nonlinear part N() and linear part C. The nonlinear part is found using Lax-Friedrichs flux. Then we apply linearisation that gives us linear operator P| u 0 := N u 0 + C near the solution u 0 . The linearisation of N() and linear matrix C becomes split in diagonal and off-diagonal parts P := L + D + U: Here λ max is the spectral radius of the Jacobi matrix that can be estimated by the maximum absolute value of the spectrum in the linearized system. Observe that matrices L, U and D consist  (1) dimension. For gas dynamics problems these are blocks of 5 × 5 and 3 × 3 for incompressible flows. First we add shift to the system by the scalar parameter σ that influences only diagonal block matrix D, obtaining modified diagonalD = D − σE.
Then for each smoother in the GM G(P, r) we perform the following sweeps: In this case Krylov residuals r n+1 are constructed using preconditioner: Sweeps (24) are performed in the menner of LU-SGS method analogous to the suggested approach in [26,27]. First we decompose our computational domain into separate subdomains for each computational device (GPU or CPU) and color these domains in such manner, that adjusted boundary elements belong to domains with different coloring, say "blue" and "green". Next, we select internal and boundary elements for these domains. First, we perform sweeps in "blue" internal cells, then we perform sweeps in half of "green" internal and all boundary cells, then we perform sweeps in "blue" boundary cells and, finaly, use the rest of the "green" cells. Inside every cell subdomain we use red-black coloring to perform GPU or OpenMP parallelization. Boundary cells are exchanged using CUDA-aware MPI. After that we perform either restriction or prolongation on GMG and perform sweeps again on different grids. It is noted in [24] that a standard preconditioner to B is not very efficient so one must apply a tuned preconditioner formed by Sherman-Morrison formula. However we were unable to find any implementation of such preconditioner for HPC architecture and this problem is the topic of another research.

Testing problems and efficiency
In this section we give results on testing problems, efficiency estimates and examples of the method application.

Matrix Market examples
First we present some test on standard matrices. Small matrix rdb450 from matrix market collection [31] is used to illustrate the application of the method. This non-symmetric matrix is originated from the discretization of reaction-diffusion brusselator model. Different stages of the method application are illustrated in figure 3. Green and red circles represent eigenvalues that are found by the application of the described inexact matrix exponent preconditioner and IRAM. We seek least stable part of the spectrum, i.e. eigenvalues with maximum real parts. Original spectrum of the matrix and desired eigenvalues are depicted in figure 3a. First we perform inexact matrix exponent transformation, result is presented in figure 3b. Observe, that the desired eigenvalues are now largest magnitude eigenvalues (LM) and are located close to the line Re(z) = 1. Then we perform shift with σ = 1.05 transforming the whole spectrum close to the line Re(z) = 0.05. Then we inverse the transformed matrix using GMRES method (in this case it converges to machine epsilon with 4 iterations) and the rest of the spectrum is inverseed analytically for illustrative purpose. The result is presented in figure 3c. We can observe that the method successfully recovered eigenvalues near the origin. However, one more complex-conjugate eigenvalue is not recovered (marked with red circles in figure 3c). In order to recover this eigenvalue we apply rotation to the matrix transformation using (17) with rotation angle 0 < ϕ < π/8 and obtain missed eigenvalues at ϕ = π/11, shown in figure 3d with the transformed spectrum. Next we consider the stress test problem of finding four eigenvalues with largest real part for non-symmetric matrix S80PI n1, sized 4025 × 4025. This matrix has two complex-conjugate desired eigenvalues λ 1,2 = 2.0 · 10 −12 ± 1.87i and λ 3,4 = 2.0 · 10 −12 ± 1.68i. We compare the suggested transformation (11) with Cayley transformation (9). For this test we chose σ = 0.01, σ 0 = −0.01, 0 < ϕ < π/3 with step dϕ = π/50 and use m = 3 with µ = 0.3. We add 10 additional Krylov vectors for filtering in the IRAM, thus Krylov subspace dimension is 14.
Results for Cayley transformation are provided in figure 4. One can see that it is impossible to recover desired spectrum with this transformation, despite good convergence of Ritz eigenvalues. The maped spectrum (middle figure) is tuning Ritz eigenvectors to the boundary. On the other hand we are able to recover desired spectrum using suggested transformation, see figure 5. One can see, that the maped spectrum of the problem is transformed in such way, that the IRAM recovers selected eigenvalues very efficiently. Red circles conrespond to ϕ = π/6.5 and cyan circle -to ϕ = π/7.5. This process takes some time, because we need to solve the eigenvalue problem for each value fo ϕ in accordance with the introduced step size dϕ.
Provided examples show that the suggested preconditioning method can be used to recover the most dangerous part of the spectrum from the linear stability point of view.

Fluid dynamics examples
Now let us consider the acceleration achieved for multiple GPU architecture and compare acceleration and efficiency with shift and inverse transformation using preconditioners (25).
Results of convergence for different methods are provided in table 1. One can see, that the choice of parameters σ, m and µ can heavily influence the method. For example for the same problem size (say 6.4E6) we can obtain convergence for the restarted GMRES method from 1100 to 3600 iterations. unfortunately, there's no a priori recipe for selecting these parameters. This we leave for future work. Second, we can see, that the BiCGStav(L) method converges more rapidly, so in practice we use this method more often.
Verification of the method on fluid dynamics applications is performed in all problems we solved, e.g. [33], [34], [35]. For example we present here the comparison of growth rates as   [35] are presented in figure 6 and wall time for compressible flow problems [34] are presented in figure 7. We can see, that only four iterations of the full IRAM are required to achieve desired tolerance. The same rapid convergence was observed in figure 5. This is achieved by adjusting regulating parameters µ and σ for the particular problem.  We noticed that the computational cost of the suggested method is O(N ) for compressible flow problem and O(N log 2 N ) for incompressible problem using spectral method. It means that the suggested method doesn't introduce additional difficulties during computation compared to the original time stepping method proportional to the problem size. Additionally one can see that the initial stagnation is related to the incomplete GPU loading.
Comparison of the preconditioned shift and inverse method with the suggested method is shown in right figure 7 and done on the compressible problem of the Kelvin-Helmholtz instability [34]. We can notice that preconditioned shift and inverse method is more efficient for small matrices, however, starting with problem size around 1 · 10 6 the presented method becomes more efficient. This can partially be explained by less efficient preconditioner for the shift  Figure 6. Convergence of full IRAM iterations and wall time required to solve the problem of finding 6 largest real eigenvalues in spectral method for incompressible flow for the residual target 1 · 10 −10 in L ∞ -norm. Figure 7. Wall time required to solve the problem of finding 6 largest real eigenvalues using IRAM for the residual target 1 · 10 −10 in L ∞ -norm and ratio of wall time for shift -inverse method using Lax-Friedrichs and GMG preconditioner. and inverse method. One can also check that the efficiency of parallel implementation for the preconditioned shift and inverse method becomes lower with the increase of parallel nodes. This is due to the increase of PCIExpress communications.
Convergence of eigenvalues to true ones in the physical problem is the function of problem discretization. In figure 8 we show the convergence of first four leading eigenvalues to the true ones in the Kelvin-Helmholtz problem as function of domain discretization and show modulo of the first unstable eigenvector.

Conclusion
In the paper we present the method of preconditioning for the linear stability eigenvalue problem. The method is aimed at the implicitly restarted Arnoldi method and can be easily added to the existing linear analysis code for fluid dynamics problems. The main purpose of the method is to be used for stability analysis of stationary solutions on cluster based on multiple GPU parallel architecture. We show that the method has analytical convergence bound and application examples confirm this bound. We also demonstrated efficient parallel application of the method in fluid dynamics problems sized up to 2 · 10 8 elements. Our future work is to tune the method and give suggestions for the selection of tuning parameters. The method implementation is available freely under GPL licence by request to the author.