---
abstract: 'The solution of eigenproblems is often a key computational bottleneck that limits the tractable system size of numerical algorithms, among them electronic structure theory in chemistry and in condensed matter physics. Large eigenproblems can easily exceed the capacity of a single compute node, thus must be solved on distributed-memory parallel computers. We here present GPU-oriented optimizations of the ELPA two-stage tridiagonalization eigensolver (ELPA2). On top of its existing cuBLAS-based GPU offloading, we add a CUDA kernel to speed up the back-transformation of eigenvectors, which can be the computationally most expensive part of the two-stage tridiagonalization algorithm. We demonstrate the performance of this GPU-accelerated eigensolver by benchmark calculations on two hybrid CPU-GPU architectures, namely a compute cluster based on Intel Xeon Gold CPUs and NVIDIA Volta GPUs, and the Summit supercomputer based on IBM POWER9 CPUs and NVIDIA Volta GPUs. Consistent with previous benchmarks on CPU-only architectures, the GPU-accelerated two-stage solver exhibits a parallel performance superior to the one-stage solver.'
address:
- 'Department of Mechanical Engineering and Materials Science, Duke University, Durham, NC 27708, USA'
- 'Molecular Sciences Software Institute, Blacksburg, VA 24060, USA'
- 'Max Planck Computing and Data Facility, Garching D-85747, Germany'
- 'Institute of Mathematics of the Czech Academy of Sciences, Žitná 25, 115 67 Prague, Czech Republic'
- NVIDIA Corporation
- 'Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37830, USA'
author:
- 'Victor Wen-zhe Yu'
- Jonathan Moussa
- Pavel Ks
- Andreas Marek
- Peter Messmer
- Mina Yoon
- Hermann Lederer
- Volker Blum
bibliography:
- 'elpa2gpu.bib'
title: 'GPU-Acceleration of the ELPA2 Distributed Eigensolver for Dense Symmetric and Hermitian Eigenproblems'
---

Eigensolver ,dense linear algebra ,parallel computing ,high-performance computing ,GPU ,CUDA

Introduction {#sec:intro}
============

Finding the eigenvalues and eigenvectors of large dense matrices is a frequent problem in computational science and engineering. For example, in Kohn-Sham density-functional theory (KS-DFT) [@dft_hohenberg_1964; @dft_kohn_1965], the many-electron problem for the Born-Oppenheimer electronic ground state is reduced to a system of single particle equations that can be discretized into a generalized eigenproblem in the following matrix form $$\label{eq:gevp}
\boldsymbol{H} \boldsymbol{C} = \boldsymbol{S} \boldsymbol{C} \boldsymbol{\Sigma} .$$

Here the Hamiltonian matrix $\boldsymbol{H}$ and the overlap matrix $\boldsymbol{S}$ are real symmetric or complex Hermitian. The matrix $\boldsymbol{C}$ and the diagonal matrix $\boldsymbol{\Sigma}$ are the eigenvectors and eigenvalues, respectively, of this eigensystem. In the framework of KS-DFT, $\boldsymbol{C}$ and $\boldsymbol{\Sigma}$ (or the information they carry, at least) are needed for the construction of $\boldsymbol{H}$. Therefore, Eq. \[eq:gevp\] is a non-linear problem and must be solved self-consistently. It is possible, and has already been implemented in various codes [@conquest_gillan_2007; @fhiaims_blum_2009; @onetep_skylaris_2005; @siesta_soler_2002], to restrict the computational cost of the construction of $\boldsymbol{H}$ to scale linearly with respect to the system size $N$ for any semi-local and hybrid exchange-correlation functional. In contrast, the solution of a dense eigenproblem (“diagonalization”) scales as $\mathcal{O}(N^3)$, quickly growing to become prohibitive as $N$ increases to large values.

Various strategies exist to facilitate or circumvent the solution of Eq. \[eq:gevp\] [@elsi_yu_2018; @elsi_yu_2019]. (1) The key idea in a conventional dense eigensolver [@numerical_press_2007; @matrix_golub_2013] is tridiagonalization, which brings the original matrix to a tridiagonal form by a series of Householder transformations. This algorithm suffers from the inefficiency of BLAS level-2 matrix-vector operations. New algorithms such as pentadiagonalization [@eigenexa_imamura_2011] and two-stage tridiagonalization [@2stage_lang_1993; @2stage_bischof_1994; @elpa_auckenthaler_2011; @elpa_marek_2014] have been developed, leading to enhanced performance over the conventional one-stage tridiagonalization approach. (2) Iterative eigensolvers [@davidson_davidson_1975; @davidson_sleijpen_1996; @iterative_payne_1992; @iterative_kresse_1996] are commonly employed by DFT codes, particularly those based on plane-wave (PW) basis functions and pseudopotentials. In that case, because of the large number of basis functions, i.e. the large dimension of the matrices in Eq. \[eq:gevp\], needed in an accurate calculation, a direct solution of Eq. \[eq:gevp\] is rather infeasible. Iterative eigensolvers are well suited to find a small fraction of low-lying eigenstates, commensurate with the needs of a PW-based code. (3) When using spatially localized basis functions, such as linear combination of atomic orbitals (LCAO), the fraction of needed eigenpairs out of the full matrix dimension can be fairly large. In this scenario, iterative solvers no longer have an advantage over direct eigensolvers. With localized basis functions, locality in the physical system can be translated to sparsity in the $\boldsymbol{H}$ and $\boldsymbol{S}$ matrices. Methods exploiting this sparsity can be formulated as $\mathcal{O}(N) \sim \mathcal{O}(N^2)$ [@linear_goedecker_1999; @linear_bowler_2012; @linear_moussa_2019; @feast_polizzi_2009; @pexsi_lin_2013; @ntpoly_dawson_2018] by circumventing the explicit solution of Eq. \[eq:gevp\]. In particular, linear scaling algorithms in a density matrix formalism have been successfully applied to simulations of one million atoms [@million_bowler_2010; @million_vandevondele_2012]. Despite the success in extreme-scale simulations, reduced scaling methods come with a computational prefactor that is much larger than that of the $\mathcal{O}(N^3)$ diagonalization method. Moreover, the applicability and optimal performance of reduced scaling methods are often limited to some certain problem types. As of today, dense eigensolvers, with their small computational prefactor and general applicability, remain the default method in most LCAO codes. Even in PW codes, the performance of a dense eigensolver is still crucial, because at some stage of an iterative solver there will be a dense eigenproblem, the size of which scales with the number of valence electrons in the system being simulated. Therefore, any improvement made to a dense eigensolver would benefit the entire electronic structure community, and the broader field of computational science in general.

The ubiquitous adoption of graphics processing units (GPU) in high-performance computing opens up new opportunities to accelerate dense eigensolvers. A GPU device consists of hundreds to thousands of parallel cores operating at a relatively low frequency. These cores are naturally suited for parallel computational tasks, such as vector and matrix operations found in dense linear algebra. On top of that, GPUs typically have a power efficiency superior to traditional CPUs, and therefore play an important role in supercomputing towards the exascale. According to the latest release of the TOP500 list at the time of writing (November 2019) [@top500], five of the top ten machines are GPU-accelerated, including Summit, the world’s fastest computer based on IBM POWER9 CPUs and NVIDIA Tesla Volta V100 GPUs.

Eigensolvers operating on hybrid CPU-GPU machines have long been available in GPU-oriented linear algebra packages such as cuSOLVER [@cusolver] and MAGMA [@magma_tomov_2010; @magma_dongarra_2014]. These packages are designed and optimized for shared-memory host architectures. They can be very fast, but the problem size they can tackle is limited by the memory capacity of a single compute node. Fully exploiting the power of GPU-accelerated supercomputers would require a distributed-memory implementation. The MPI-parallel, distributed-memory ELPA library implements the conventional one-stage diagonalization method and the two-stage diagonalization proposed in Refs. [@2stage_lang_1993; @2stage_bischof_1994], known as the “ELPA1” and “ELPA2” solvers, respectively. Both of them have been ported to GPUs by substituting BLAS calls with the corresponding cuBLAS functions [@elpa_kus_2019a], making ELPA the only publicly available, distributed-memory, GPU-accelerated eigensolver to our knowledge. The CPU-only version of ELPA1 and ELPA2 and the GPU-accelerated version of ELPA1 and ELPA2 are hereafter referred to as CPU-ELPA1, CPU-ELPA2, GPU-ELPA1, and GPU-ELPA2, respectively.

In single-node tests published by Ks et al. [@elpa_kus_2019a; @elpa_kus_2019b], the performance of GPU-ELPA1 is better than that of GPU-ELPA2. When using 2 IBM POWER8 CPUs (24 cores in total) and 4 NVIDIA Pascal P100 GPUs, GPU-ELPA1 delivers up to 11.9x performance boost compared to CPU-ELPA1 using 24 CPU cores. The efficiency of CPU-ELPA2 on various distributed-memory CPU platforms [@elpa_marek_2014; @elpa_gutheil_2014; @elpa_cook_2018] indicates that GPU-ELPA2 may outperform GPU-ELPA1 on multiple nodes. Historically, the GPU porting of ELPA2 dates back to 2013, when Peter Messmer from NVIDIA programmed the first version of GPU-ELPA2. Then the code was refactored and merged into the mainline version of ELPA, and has been available in released versions of the ELPA eigensolver library since 2016. In this paper, we report our latest optimizations and developments that enable a performance improvement on distributed-memory, GPU-accelerated architectures. Specifically, several synchronizations and memory transfers between CPUs and GPUs have been optimized. Additionally, some kernels in one of the major computational steps, the tridiagonal-to-banded back-transformation of eigenvectors (Eq. \[eq:2stage\_bkwd1\]), have been rewritten. The rest of this paper is organized as follows. First, we briefly review the two-stage diagonalization algorithm, in particular the tridiagonal-to-banded back-transformation of eigenvectors and its CPU implementation in the ELPA library. Next, we outline the GPU acceleration strategies employed in GPU-ELPA2, and elaborate on our CUDA implementation of the tridiagonal-to-banded back-transformation of eigenvectors, which is essentially a GPU extension of the algorithm in Refs. [@elpa_auckenthaler_2011; @elpa_auckenthaler_2013]. Finally, we demonstrate the performance and scalability of the GPU-ELPA2 solver by a set of benchmarks on two GPU-accelerated computers, namely the Talos cluster at Max Planck Computing and Data Facility in Garching, Germany, based on Intel Xeon Gold CPUs and NVIDIA Volta GPUs, and the Summit supercomputer at Oak Ridge National Laboratory in Tennessee, USA, based on IBM POWER9 CPUs and NVIDIA Volta GPUs.

Two-Stage Tridiagonalization in ELPA2 {#sec:review_2stage}
=====================================

Overview of the Two-Stage Tridiagonalization {#subsec:cpu_2stage}
--------------------------------------------

The textbook procedure [@numerical_press_2007; @matrix_golub_2013] to solve a dense generalized eigenproblem, like the one in Eq. \[eq:gevp\], first computes the Cholesky factorization of $\boldsymbol{S}$ $$\label{eq:cholesky}
\boldsymbol{S} = \boldsymbol{L} \boldsymbol{L}^* ,$$

then uses $\boldsymbol{L}$ to transform Eq. \[eq:gevp\] to a standard eigenproblem $$\label{eq:evp}
\boldsymbol{\tilde{H}} \boldsymbol{\tilde{C}} = \boldsymbol{\tilde{C}} \boldsymbol{\Sigma} .$$

$\boldsymbol{\tilde{H}}$ is $$\label{eq:gevp2evp}
\boldsymbol{\tilde{H}} = \boldsymbol{L}^{-1} \boldsymbol{H} (\boldsymbol{L}^*)^{-1} ,$$

and the eigenvectors $\boldsymbol{\tilde{C}}$ must be back-transformed in order to retrieve the eigenvectors of Eq. \[eq:gevp\], i.e. $$\label{eq:evp2gevp}
\boldsymbol{C} = (\boldsymbol{L}^*)^{-1} \boldsymbol{\tilde{C}} .$$

The direct solution of Eq. \[eq:evp\] is based on tridiagonalization, that is, the full matrix $\boldsymbol{\tilde{H}}$ is transformed to a tridiagonal matrix $\boldsymbol{T}$. This is typically accomplished by individual Householder transformations, which take the shape of matrix-vector operations. Eigenvalues and eigenvectors of $\boldsymbol{T}$ can be easily (compared to the original problem) solved. Then, eigenvectors of $\boldsymbol{T}$ are back-transformed to obtain eigenvectors of $\boldsymbol{\tilde{H}}$. This algorithm is adopted by a variety of dense linear algebra packages, such as LAPACK [@lapack_anderson_1999] targeting sequential and shared-memory parallel architectures; cuSOLVER [@cusolver] and MAGMA [@magma_tomov_2010; @magma_dongarra_2014] targeting shared-memory architectures with GPU accelerators; ScaLAPACK [@scalapack_blackford_1997] and Elemental [@elemental_poulson_2013] targeting distributed-memory parallel architectures. However, as mentioned above, the tridiagonalization step makes extensive use of memory-bound, BLAS level-2 matrix-vector operations, whose performance is limited on modern computer architectures.

The two-stage tridiagonalization algorithm proposed by Bischof, Sun, and Lang [@2stage_lang_1993; @2stage_bischof_1994] is an established alternative to the conventional one-stage method. As shown in Eq. \[eq:2stage\] below and further illustrated in Fig. \[fig:2stage\], the tridiagonalization of the full matrix $\boldsymbol{\tilde{H}}$ is carried out in two transformations. The first transformation $\boldsymbol{P}$ reduces $\boldsymbol{\tilde{H}}$ to a banded matrix $\boldsymbol{B}$, and the second transformation $\boldsymbol{Q}$ reduces $\boldsymbol{B}$ to a tridiagonal matrix $\boldsymbol{T}$. The eigenvalues $\boldsymbol{\Sigma}$ and eigenvectors $\boldsymbol{X}$ of $\boldsymbol{T}$ are solved as done in the one-stage method. The back-transformation of eigenvectors is also carried out in two steps. $\boldsymbol{X}$ is first back-transformed to $\boldsymbol{Y}$, the eigenvectors of $\boldsymbol{B}$, then to $\boldsymbol{\tilde{C}}$, the eigenvectors of $\boldsymbol{\tilde{H}}$.

\[eq:2stage\] $$\begin{aligned}
\boldsymbol{B} & = \boldsymbol{P} \boldsymbol{\tilde{H}} \boldsymbol{P}^* , \label{eq:2stage_fwd1} \\
\boldsymbol{T} & = \boldsymbol{Q} \boldsymbol{B} \boldsymbol{Q}^* , \label{eq:2stage_fwd2} \\
\boldsymbol{T} \boldsymbol{X} & = \boldsymbol{X} \boldsymbol{\Sigma} , \label{eq:2stage_solve} \\
\boldsymbol{Y} & = \boldsymbol{Q}^* \boldsymbol{X} , \label{eq:2stage_bkwd1} \\
\boldsymbol{\tilde{C}} & = \boldsymbol{P}^* \boldsymbol{Y} . \label{eq:2stage_bkwd2}\end{aligned}$$

![image](fig/2stage.png){width="99.00000%"}

This two-stage tridiagonalization approach is implemented in several linear algebra software packages [@elpa_auckenthaler_2011; @elpa_marek_2014; @magma_tomov_2010; @lapack_anderson_1999; @mkl_arturov_2018], including a high-performance distributed-memory implementation in the ELPA library [@elpa_auckenthaler_2011; @elpa_marek_2014]. The introduction of the banded matrix stage leads to faster computation compared to the one-stage method, for the transformation in Eq. \[eq:2stage\_fwd1\] mostly involves highly efficient BLAS level-3 matrix-matrix operations, and the transformation in Eq. \[eq:2stage\_fwd2\] only works on a sparse banded matrix $\boldsymbol{B}$ instead of a full matrix. The solution of Eq. \[eq:2stage\_solve\] is accelerated in ELPA by extending the divide-and-conquer symmetric tridiagonal eigensolver [@dc_cuppen_1980; @dc_gu_1995; @dc_tisseur_1999] such that unwanted eigenvectors are not computed [@elpa_auckenthaler_2011; @elpa_auckenthaler_2013]. Regarding the back-transformation of eigenvectors, it is rather difficult to directly use BLAS routines for Eq. \[eq:2stage\_bkwd1\]. Manually optimized “kernels” written in architecture-specific instruction sets are employed for this particular step [@elpa_auckenthaler_2011; @elpa_kus_2019b; @elpa_auckenthaler_2013]. The ELPA2 solver is highly scalable on massively parallel, distributed-memory architectures. It avoids global MPI communications by using a 2-dimensional (2D) process grid and restricting the communication to take place within either the row direction or the column direction. Depending on the size of the matrix, ELPA2 scales to at least tens of thousands of CPU cores [@elpa_marek_2014; @elpa_gutheil_2014; @elpa_cook_2018].

Tridiagonal-to-Banded Back-Transformation of Eigenvectors {#subsec:cpu_step4}
---------------------------------------------------------

We now summarize the tridiagonal-to-banded back-transformation algorithm and its CPU implementation in ELPA2. The reader is also referred to Refs. [@2stage_lang_1993; @2stage_bischof_1994; @elpa_auckenthaler_2011; @elpa_auckenthaler_2013]. ELPA2 relies on the “bulge chasing” algorithm [@2stage_lang_1993; @2stage_bischof_1994] to reduce a banded matrix $\boldsymbol{B}$ to a tridiagonal matrix $\boldsymbol{T}$. Let $N$ and $b$ denote the dimension and semi-bandwidth of $\boldsymbol{B}$, respectively. The banded-to-tridiagonal transformation is done in $N-2$ stages, with $(N-i)/b$ sweeps in the $i^\text{th}$ stage. The first sweep of the $i^\text{th}$ stage reduces the $i^\text{th}$ column of $\boldsymbol{B}$ to the target tridiagonal form, at the same time introducing fill-ins (“bulges”) to the remainder of $\boldsymbol{B}$. From the second sweep on, the first column of the fill-ins introduced in the previous sweep is eliminated, while introducing new fill-ins further down the matrix. The transformation applied in the $j^\text{th}$ sweep of the $i^\text{th}$ stage can be written as $$\label{eq:householder}
\boldsymbol{Q}_{(i,j)} = \boldsymbol{Q}^*_{(i,j)} = \boldsymbol{I} - \tau_{(i,j)} \boldsymbol{v}_{(i,j)} \boldsymbol{v}_{(i,j)}^* ,$$

where $\boldsymbol{Q}_{(i,j)}$ is a Householder transformation matrix; $\boldsymbol{I}$ is the identity matrix; the scalar $\tau_{(i,j)}$ and vector $\boldsymbol{v}_{(i,j)}$ are computed following the standard Householder method [@numerical_press_2007; @householder_householder_1958]. Most Householder vectors $\boldsymbol{v}_{(i,j)}$ have a length equal to $b$, except that $\boldsymbol{v}_{(i,j)}$ generated in the last sweep of each stage may be shorter. Matrix-vector operations are still needed in order to apply these Householder transformations, but the computational cost is much less than in the one-stage tridiagonalization algorithm, as the vectors are much shorter. The left panel of Fig. \[fig:step4\] shows where vectors $\boldsymbol{v}_{(i,j)}$ are generated with an example of $N=17$, $b=4$, and $N-2=15$ stages of bulge chasing. In the actual code, the Householder matrices $\boldsymbol{Q}_{(i,j)}$ are never explicitly constructed. Instead, $\tau_{(i,j)}$ and $\boldsymbol{v}_{(i,j)}$ are stored and used for the transformations in Eqs. \[eq:2stage\_fwd2\] and \[eq:2stage\_bkwd1\], where $\boldsymbol{Q}$ is the product of all $\boldsymbol{Q}_{(i,j)}$.

![image](fig/step4.png){width="99.00000%"}

After the tridiagonal eigenproblem (Eq. \[eq:2stage\_solve\]) is solved, we get the eigenvalues $\boldsymbol{\Sigma}$ and eigenvectors $\boldsymbol{X}$ of the tridiagonal matrix $\boldsymbol{T}$. The main task of the tridiagonal-to-banded back-transformation is to apply all $\boldsymbol{Q}^*_{(i,j)}$ to $\boldsymbol{X}$ (Eq. \[eq:2stage\_bkwd1\]). It is obvious that the eigenvectors can be back-transformed independently, leading to a trivial parallelism over eigenvectors. Moreover, as the Householder vectors $\boldsymbol{v}_{(i,j)}$ are shorter than the eigenvectors $\boldsymbol{X}$, a single transformation $\boldsymbol{Q}^*_{(i,j)}$ only alters a few rows of $\boldsymbol{X}$ where $\boldsymbol{v}_{(i,j)}$ has non-zero values. For instance, the Householder vector $\boldsymbol{v}_{(1,4)}$ at the left bottom corner of Fig. \[fig:step4\] only alters the last four rows of $\boldsymbol{X}$. This leads to another level of parallelism to be exploited within an individual eigenvector.

In ELPA2, the $N_p$ processes are organized in an $N_{pr}$ by $N_{pc}$ grid. The $N_{ev}$ eigenvectors are uniformly distributed across the $N_{pc}$ process columns. Within a process column, the eigenvectors are distributed across the $N_{pr}$ processes in a block manner. Each process applies a series of Householder transformations to its local part of the eigenvector matrix $\boldsymbol{X}$. The colored part in Fig. \[fig:step4\] shows an example of the local computation by a process. This process applies four Householder transformations, $\boldsymbol{Q}^*_{(8,1)}$, $\boldsymbol{Q}^*_{(7,1)}$, $\boldsymbol{Q}^*_{(6,1)}$, and $\boldsymbol{Q}^*_{(5,1)}$, to its local part of the eigenvector matrix, referred to as $\boldsymbol{X}_\text{local}$ hereafter. Apparently, $\boldsymbol{Q}^*_{(4,1)}$, $\boldsymbol{Q}^*_{(3,1)}$, and $\boldsymbol{Q}^*_{(2,1)}$ also modify the top part of $\boldsymbol{X}_\text{local}$. These three transformations are however applied by the upper adjacent process in the same process column. Therefore, the top three rows of $\boldsymbol{X}_\text{local}$ must be exchanged with the upper adjacent process. Likewise, the bottom three rows of $\boldsymbol{X}_\text{local}$ must be exchanged with the lower adjacent process in the same process column. The middle part of $\boldsymbol{X}_\text{local}$ is not involved in any data exchange. It appears in Fig. \[fig:step4\] that almost the entire $\boldsymbol{X}_\text{local}$ matrix needs to be exchanged between adjacent process rows, as there is only one row in the middle part. However, the height of the middle part actually increases with the matrix size $N$, whereas the height of the top and bottom parts can never exceed the semi-bandwidth $b$. Therefore, the amount of data that needs to be exchanged, i.e., data in the top and bottom parts, is limited by $b$, usually accounting for a small fraction of $\boldsymbol{X}_\text{local}$.

GPU Acceleration of ELPA2 {#sec:gpu_2stage}
=========================

In ELPA2, the two-stage tridiagonalization is implemented as five separate subroutines corresponding to the five steps in Eq. \[eq:2stage\]. The input, output, and internal working matrices of ELPA2 are distributed across CPU cores. ELPA2 relies on its own explicit MPI calls to handle distributed linear algebra operations. The efficient MPI communication pattern in the CPU version of ELPA2 is not altered in the GPU version, where GPU acceleration is mainly realized by substituting local BLAS calls with the corresponding cuBLAS functions, as is done for the full-to-banded transformation (Eq. \[eq:2stage\_fwd1\]), the solution of the tridiagonal system (Eq. \[eq:2stage\_solve\]), and the banded-to-full back-transformation (Eq. \[eq:2stage\_bkwd2\]). The tridiagonal-to-banded back-transformation (Eq. \[eq:2stage\_bkwd1\]) is GPU accelerated by a CUDA implementation of Eq. \[eq:householder\]. The banded-to-tridiagonal transformation (Eq. \[eq:2stage\_fwd2\]) has not been ported to GPUs, because of its low computational cost as shown in Fig. \[fig:decomp5\_summit\].

GPU Offloading via cuBLAS {#subsec:gpu_cublas}
-------------------------

The API of cuBLAS is designed to be almost identical to that of the standard CPU BLAS, making cuBLAS-based GPU offloading straightforward. Here we only comment on two technical aspects. First, the communication and synchronization between CPUs and GPUs should be avoided as much as possible. Before calling a cuBLAS function, the input arrays must reside on the GPU memory, which often requires a copy of the data from CPU to GPU. In order for the CPU to access the result of cuBLAS, another copy from GPU to CPU is needed. In our GPU porting of ELPA2, CPU-GPU memory transfers are reduced to minimum by leaving data on the GPU as long as possible. Most often, GPU data is copied to the CPU in order to participate in an MPI communication call. In the version discussed in this paper, we have not yet explored GPU-aware MPI to directly communicate data on the GPU. In each of the GPU-accelerated subroutines, the allocation of GPU memory, which implies a CPU-GPU synchronization, is performed before the main work begins by precomputing the size of the allocation. The allocated GPU memory is reused wherever possible, and is deallocated after the main work finishes. Avoiding frequent allocation and deallocation of GPU memory in a loop or in a recursive routine greatly reduces the amount of implicit CPU-GPU synchronization.

Second, the CPU code of the GPU-accelerated version of ELPA2 operates in a pure MPI mode without threading, i.e., one MPI task for each CPU core. As most (if not all) mainstream computers today have more CPU cores than GPU devices, several MPI tasks would have to share one GPU device. When the size of an eigenproblem is relatively small, the amount of work assigned to each individual MPI task may not be able to fully saturate the GPU. In such cases, the NVIDIA Multi-Process Service (MPS) transparently allows work from different MPI tasks to be executed concurrently on the GPU, increasing the overall GPU utilization. It is thus recommended to use ELPA with MPS switched on.

CUDA Implementation of Parallel Householder Transformations {#subsec:gpu_step4}
-----------------------------------------------------------

In this section, we present our CUDA implementation of the local computation in Fig. \[fig:step4\], which is the key step in the tridiagonal-to-banded back-transformation (Eq. \[eq:2stage\_bkwd1\]). In a CUDA program, the large number of GPU cores are arranged into a grid of blocks, each of which in turn comprises a grid of threads. All the GPU cores work in a single instruction, multiple thread (SIMT) fashion, i.e., a single instruction is simultaneously executed by multiple threads with different data [@cuda_nickolls_2008]. As explained in Sec. \[subsec:cpu\_step4\], in a given step of the tridiagonal-to-banded back-transformation, a process is responsible for applying $N_v$ Householder transformations to its local eigenvector matrix $\boldsymbol{X}_\text{local}$, which has $N_R$ rows and $N_C$ columns. A Householder transformation defined by $\tau$ and $\boldsymbol{v}$ is applied to an eigenvector $\boldsymbol{x}$ by $$\label{eq:kernel}
(\boldsymbol{I} - \tau \boldsymbol{v} \boldsymbol{v}^*) \boldsymbol{x} = \boldsymbol{x} - \tau \boldsymbol{v} (\boldsymbol{v}^* \boldsymbol{x}) .$$

The $N_v$ Householder transformations are applied sequentially. The order in which they are applied in the tridiagonal-to-banded back-transformation is the reverse of the order in which they are generated in the banded-to-tridiagonal transformation. Consider the example in Fig. \[fig:step4\] again, the first transformation would be $\boldsymbol{Q}^*_{(8,1)}$, then $\boldsymbol{Q}^*_{(7,1)}$, $\boldsymbol{Q}^*_{(6,1)}$, and finally $\boldsymbol{Q}^*_{(5,1)}$. It is obvious that from one transformation to the next, the rows of $\boldsymbol{X}_\text{local}$ modified by the transformation are shifted upward by one row. This is also seen in Fig. \[fig:step4\], where $\boldsymbol{Q}^*_{(8,1)}$ modifies the 4$^\text{th}$ to 7$^\text{th}$ rows of $\boldsymbol{X}_\text{local}$, $\boldsymbol{Q}^*_{(7,1)}$ modifies the 3$^\text{rd}$ to 6$^\text{th}$ rows of $\boldsymbol{X}_\text{local}$, and so forth. Each individual $\boldsymbol{Q}^*_{(i,j)}$ actually only modifies $b$ rows of $\boldsymbol{X}_\text{local}$.

In order to map this local computation to the GPU cores, we choose a 1D block grid with a 1D thread grid within each block. The CUDA kernel is launched with $N_C$ blocks and $b$ threads, so that each block works on an eigenvector, and each thread works on an element of this eigenvector. Specifically, the Householder transformations are applied to $\boldsymbol{X}_\text{local}$ as follows:

1.  Copy $\boldsymbol{X}_\text{local}$, as well as all $\boldsymbol{v}$ and $\tau$ that ever update $\boldsymbol{X}_\text{local}$, from CPU to GPU. Each GPU block is responsible for one column of $\boldsymbol{X}_\text{local}$, denoted as $\boldsymbol{x}$.

2.  Compute Eq. \[eq:kernel\] for the first Householder transformation.

    1.  Compute dot product $\boldsymbol{v}^* \boldsymbol{x}$. Only $b$ elements of $\boldsymbol{x}$ contribute to the dot product. Each of the $b$ GPU threads loads one element of $\boldsymbol{v}$ and one element of $\boldsymbol{x}$, multiplies them together, and stores the thread-local result in shared memory. The final dot product is obtained by a parallel reduction involving all the threads in a block.

    2.  Update $\boldsymbol{x}$ by $\boldsymbol{x} = \boldsymbol{x} - \tau (\boldsymbol{v}^* \boldsymbol{x}) \boldsymbol{v}$. Again, only $b$ elements of $\boldsymbol{x}$ are modified. Since the dot product $\boldsymbol{v}^* \boldsymbol{x}$ has already been computed and $\tau$ is a scalar, this update can be done in a straightforward element-wise fashion, i.e., each thread updates one element of $\boldsymbol{x}$.

3.  Among the $b$ elements of $\boldsymbol{x}$ that are updated by the first Householder transformation, only the lowest element will not be affected by the next Householder transformation. Before applying the next transformation, the last thread in each block writes its element of $\boldsymbol{x}$ into $\boldsymbol{X}_\text{local}$. Then thread $i$ ($i>=1$) takes the element of $\boldsymbol{x}$ from thread ($i-1$), while thread 0 loads a new element from $\boldsymbol{X}_\text{local}$.

4.  Now all threads are ready to repeat steps 1, 2, and 3 for the next Householder transformation.

5.  The CUDA kernel finishes when all transformations are applied. The updated $\boldsymbol{X}_\text{local}$ is then copied back to CPU.

An example of this procedure is given in Fig. \[fig:kernel\], where four Householder transformations are applied to a local eigenvector matrix with $N_R=7$ rows and $N_C=6$ columns. The CUDA kernel is launched with six blocks and four threads, as indicated by the block id and thread id (both zero-based) in the figure. Throughout the execution of the kernel, block $i$ is responsible for the $(i+1)^\text{th}$ eigenvector. The four Householder transformations need to be applied in reverse order, i.e. from right to left. In the application of the rightmost transformation, threads 0, 1, 2, 3 in each block are responsible for the 4$^\text{th}$, 5$^\text{th}$, 6$^\text{th}$, and 7$^\text{th}$ elements, respectively, of the eigenvector this block is responsible for. After computing Eq. \[eq:kernel\], thread 3 immediately writes the 7$^\text{th}$ element back to the $\boldsymbol{X}_\text{local}$ array. Threads 1, 2, 3 take the updated 4$^\text{th}$, 5$^\text{th}$, and 6$^\text{th}$ elements from threads 0, 1, 2, respectively. Then thread 0 loads the 3$^\text{rd}$ element of the eigenvector from $\boldsymbol{X}_\text{local}$, and all the threads are ready for the next transformation. It is obvious that all the transformations can be applied in exactly the same way. After the kernel finishes, the final $\boldsymbol{X}_\text{local}$ is copied back from GPU to CPU.

![image](fig/kernel.png){width="99.00000%"}

Performance and Scalability {#sec:results}
===========================

The performance of the GPU-accelerated ELPA2 solver is benchmarked on the Talos cluster at Max Planck Computing and Data Facility and the Summit supercomputer at Oak Ridge National Laboratory. Each node of Talos has 2 Intel Xeon Gold 6148 CPUs (40 cores in total) and 2 NVIDIA Tesla Volta V100 (each has 32 GB high-bandwidth memory, double precision peak 7.0 TFLOP/s, PCIe 32 GB/s interconnect) GPUs. Benchmarks presented in Fig. \[fig:all\_talos\] are performed on Talos. The ELPA code is compiled with the Intel 19.0.5 compiler suite, Intel MPI 2019.5, Intel MKL 2019.5, and CUDA 10.1. Each node of Summit has 2 IBM POWER9 CPUs (44 cores in total, of which 42 are for running applications) and 6 NVIDIA Tesla Volta V100 (each has 16 GB high-bandwidth memory, double precision peak 7.8 TFLOP/s, NVLink 300 GB/s interconnect) GPUs. Benchmarks presented in Figs. \[fig:all\_summit\], \[fig:decomp3\_summit\], \[fig:decomp5\_summit\] are performed on Summit. The ELPA code is compiled with the IBM XL 16.1.1 compiler suite, IBM Spectrum MPI 10.3, IBM ESSL 6.1, and CUDA 10.1.

Fig. \[fig:all\_talos\] shows the total time to solution of CPU-ELPA1, CPU-ELPA2, GPU-ELPA1, and GPU-ELPA2 on the Talos cluster. All eigenvalues and eigenvectors of randomly generated real and complex matrices of size $N$ = 40,000 to 100,000 are computed with up to 64 nodes. As already demonstrated in published benchmarks [@elpa_marek_2014; @elpa_gutheil_2014; @elpa_cook_2018; @elpa_auckenthaler_2013], CPU-ELPA2 greatly outperforms CPU-ELPA1 in terms of performance and scalability. The performance difference between GPU-ELPA1 and GPU-ELPA2 is rather small. For small node counts, GPU-ELPA1 is marginally faster than GPU-ELPA2, which agrees with the published single-node tests [@elpa_kus_2019b]. What has not been tested previously is the performance on multiple nodes. It turns out that GPU-ELPA2 becomes faster than GPU-ELPA1 as the node count increases. The crossover point depends on the dimension of the problem, e.g., 4 nodes in the $N$ = 40,000 real case, and 16 nodes in the $N$ = 100,000 complex case. When using 64 nodes, the speedup of GPU-ELPA2 over GPU-ELPA1, averaged over all matrix sizes in Fig. \[fig:all\_talos\], is 2.2x.

In Fig. \[fig:all\_talos\], the speedup enabled by the GPUs ranges from no speedup at all to 3.3x. Three general trends emerge: (1) For the same matrix size, the speedup becomes smaller as the node count increases. For small $N$, CPU-ELPA2 can even be faster than the GPU-accelerated solvers, thanks to the near-optimal strong scaling of CPU-ELPA2. (2) For the same number of nodes, the speedup becomes larger as the matrix size increases. (3) For the same node count and the same matrix size, the speedup is larger for complex matrices than for real matrices.

![image](fig/all_talos.png){width="99.00000%"}

The timing experiment in Fig. \[fig:all\_talos\] is repeated on the Summit supercomputer. The results are shown in Fig. \[fig:all\_summit\]. CPU-ELPA1 is omitted for simplicity. On Summit, the GPU-accelerated solvers GPU-ELPA1 and GPU-ELPA2 are always faster than the CPU-ELPA2 solver, with a maximum speedup of over 20x. The speedup of GPU-ELPA2 over CPU-ELPA2 remains 2.2x ($N$ = 40,000 real) to 6.7x ($N$ = 100,000 complex) even for 64 nodes. For the same matrix size and the same node count, the speedup on Summit appears greater than on Talos, which can be partially attributed to the difference in hardware. Summit has 6 GPUs per node, whereas Talos only has two GPUs per node. Data transfers on Summit take advantage of the NVLink technology [@nvlink_foley_2017] for high-bandwidth interconnect between CPUs and GPUs. Besides, a high-performance CPU kernel for the tridiagonal-to-banded back-transformation, written in AVX512 instructions [@elpa_kus_2019b], is employed for the Intel Xeon Gold CPUs on Talos, rendering better performance of CPU-ELPA2 on Talos. Overall, we therefore observe that the absolute per-node timings of CPU-ELPA2 in Figs. \[fig:all\_talos\] and \[fig:all\_summit\] are already lower on Talos than on Summit. This difference in the CPU-only results, which are the baseline of the reported speedups, probably exaggerates the comparison of GPU-ELPA2 and CPU-ELPA2 on Summit somewhat, relative to Talos. Nevertheless, the three trends summarized from the tests on Talos are still valid on Summit, that is, the GPU speedup is larger for (1) larger matrix size, (2) fewer nodes, and (3) solving a complex problem instead of a real one. On both computers, the strong scaling of the GPU solvers is never as good as that of the CPU solvers. This can be explained by the workload assigned to the individual nodes. When using a large number of nodes or solving a small matrix, the workload on each node becomes so little that the many GPUs cannot be saturated, and the cost of CPU-GPU communications cannot be amortized. In contrast, when solving a large matrix or using a small number of nodes, a large amount of local work is offloaded to the GPUs, resulting in a significant speedup.

![image](fig/all_summit.png){width="99.00000%"}

Note that in Figs. \[fig:all\_talos\] and \[fig:all\_summit\], all eigenvalues and eigenvectors are computed. In applications such as materials simulations based on KS-DFT, only a portion of the eigenspectrum, e.g., typically 20% to 60% for LCAO basis sets, is of interest. In this case, the advantage of ELPA2 over ELPA1 should be more significant, as the computational complexity of the back-transformation is proportional to the number of eigenvectors to be calculated. This is demonstrated in Fig. \[fig:decomp3\_summit\], where the total timings (red circles) of the GPU-ELPA1 (solid) and GPU-ELPA2 (dashed) are decomposed into the forward tridiagonalization (blue squares), the solution of the tridiagonal problem (yellow diamonds), and the backward transformation (violet triangles). Two test cases, namely $N$ = 40,000 real and $N$ = 100,000 real, are shown as examples for a single node on Summit, i.e., the worst-case scenario for ELPA2. The two-stage tridiagonalization in ELPA2 is always faster than the one-stage tridiagonalization in ELPA1 by a factor of approximately three. The back-transformation accounts for a small fraction of the total time in ELPA1, but is the most time-consuming part in ELPA2 when all the eigenvectors are computed. When computing fewer eigenvectors, the burden of the two-stage back-transformation in ELPA2 can be greatly alleviated, making ELPA2 more favorable than ELPA1.

![image](fig/decomp3_summit.png){width="99.00000%"}

Given that the optimal performance may be achieved with GPU-ELPA1, GPU-ELPA2, or CPU-ELPA2, depending on the specifics of the problem and the architecture, we highlight the auto-tuning feature in the ELPA library [@elpa_kus_2019b]. When ELPA is called repeatedly, like e.g. in a self-consistent KS-DFT calculation, this auto-tuning feature automatically iterates over possible combinations of solvers and runtime parameters. The best solver for a given problem can be identified and utilized without any additional input from the user.

In Fig. \[fig:decomp5\_summit\], we further decompose the timings of CPU-ELPA2 and GPU-ELPA2 into the five computational steps defined in Eq. \[eq:2stage\]. Again, the $N$ = 40,000 real and $N$ = 100,000 real tests are shown as examples. Steps that have been GPU-accelerated display an excellent speedup, namely 5.8x, 6.1x, 17.1x, and 11.3x (averaged over all data points in Fig. \[fig:all\_summit\]) for the full-to-banded transformation (Eq. \[eq:2stage\_fwd1\]), the solution of the tridiagonal problem (Eq. \[eq:2stage\_solve\]), the tridiagonal-to-banded back-transformation (Eq. \[eq:2stage\_bkwd1\]), and the banded-to-full back-transformation (Eq. \[eq:2stage\_bkwd2\]), respectively. The banded-to-tridiagonal transformation step (Eq. \[eq:2stage\_fwd2\]) is not yet GPU-accelerated, as it never stands as a bottleneck. The tridiagonal-to-banded back-transformation, which uses the newly developed CUDA kernel described in Sec. \[subsec:gpu\_step4\], shows a strong scaling that is close to ideal. The scaling of the full-to-banded transformation and the banded-to-full back-transformation is not as good as the other steps. These two steps limit the overall parallel efficiency of GPU-ELPA2, therefore they would be the first target for further algorithmic and technical optimization.

![image](fig/decomp5_summit.png){width="99.00000%"}

Conclusions {#sec:conclude}
===========

In this paper, we report our GPU-oriented optimizations of the two-stage tridiagonalization eigensolver ELPA2 in the ELPA library, the only publicly available, distributed-memory, GPU-accelerated eigensolver library at the time of writing. The local BLAS operations in ELPA2 are offloaded to GPUs via the cuBLAS library. The tridiagonal-to-banded back-transformation of eigenvectors, which cannot be easily written as BLAS operations, is GPU-accelerated by a CUDA implementation. The overall performance of the GPU-accelerated ELPA2 solver is promising. It delivers a significant performance boost over the CPU-only version of ELPA2, as demonstrated by benchmarks on the Talos and Summit computers. Owing to the advanced two-stage tridiagonalization algorithm, the parallel scaling of the GPU-ELPA2 solver is superior to that of the GPU-ELPA1 solver. Based on an analysis of the individual computational steps in ELPA2, we identify the full-to-banded transformation and banded-to-full back-transformation steps as the next target for future optimization, as their strong scaling is not yet optimal. Nevertheless, the GPU-ELPA2 solver in its current form already unlocks a significant potential of GPU computations.

Acknowledgments {#sec:thanks .unnumbered}
===============

This research was supported by the National Science Foundation (NSF) under Award 1450280. Yu was supported by a fellowship from the Molecular Sciences Software Institute under NSF Award 1547580. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC05-00OR22725. We thank Dr. Christoph Angerer and Scott Biersdorff from NVIDIA, Dr. Björn Lange, Dr. William Huhn, and Dr. Wenhui Mi from Duke University, Dr. Jack Deslippe, Dr. Mauro Del Ben, and Dr. Charlene Yang from National Energy Research Scientific Computing, Dr. Eduardo D’Azevedo from Oak Ridge Leadership Computing Facility (OLCF), Dr. Markus Rampp from Max Planck Computing and Data Facility, Garching, Germany, and Dr. Christian Carbogno and Dr. Matthias Scheffler from Fritz Haber Institute of the Max Planck Society, Berlin, Germany for fruitful discussions. Part of the optimization presented in this paper was the outcome of the 2019 OLCF GPU Hackathon. We gratefully acknowledge the organizer Dr. Thomas Papatheodore (OLCF) and mentors Brent Leback (NVIDIA) and Dr. Markus Eisenbach (OLCF) for their advice during this event.
