ESSEX: Equipping Sparse Solvers For Exascale

,

were guided by quantum physics applications, especially from the field of topological insulator-or graphene-based systems.For these, ScaMaC, a scalable matrix generation framework for a broad set of quantum physics problems, was developed.As the central software core of ESSEX, the PHIST library for sparse systems of linear equations and eigenvalue problems has been established.It abstracts algorithmic developments from low-level optimization.Finally, central ESSEX software components and solvers have demonstrated scalability and hardware efficiency on up to 256 K cores using million-way process/thread-level parallelism.

Introduction
The efficient solution of linear systems or eigenvalue problems involving large sparse matrices has been an active research field in parallel and high performance computing for many decades.Software packages like Trilinos [33] or PETSc [9] have been developed to great maturity, and algorithmic improvements were accompanied by advances in programming abstractions addressing, e.g., node-level heterogeneity (cf.Kokkos [19]).Completely new developments such as Ginkgo1 are rare and do not focus on large-scale applications or node-level efficiency.
Despite projections from the late 2000s, hardware architectures have not developed away from traditional clustered multicore systems.However, a clear trend of increased node-level parallelism and heterogeneity has been observed.Although several new architectures entered the field (and some vanished again), the basic concepts of core-level code execution and data parallelism have not changed.This is why the MPI+X concept is still a viable response to the challenge of hardware diversity.
Performance analysis of highly parallel code typically concentrated on scalability, but provably optimal node-level performance was rarely an issue.Moreover, strong abstraction boundaries between linear algebra building blocks, solvers, and applications made it hard to get a holistic view on a minimization of time to solution, encompassing optimizations in the algorithmic and implementation dimensions.
In this setting, the ESSEX project took the opportunity to start from a clean slate, deliberately breaking said abstraction boundaries to investigate performance bottlenecks together with algorithmic improvements from the core to the highly parallel level.Driven by the targeted application fields, bespoke solutions were developed for selected algorithms and applications.The experience gained in the development process will lead the way towards more generic approaches rather than compete with established libraries in terms of generality.The overarching motif was a consistent performance engineering process that coordinated all performancerelevant activities across the different software layers [1, 3, 4, 20, 52-54, 56, 62, 64].
Consequently, the ESSEX parallel building blocks layer implemented in the GHOST library [55] supports MPI+X, with X being a combination of nodelevel programming models able to fully exploit hardware heterogeneity, functional parallelism, and data parallelism.Despite fluctuations in hardware architectures and new programming models hitting the market every year, OpenMP or CUDA is still the most promising and probably most sustainable choice for X, and ESSEX-II adhered to it.In addition, engineering highly specialized kernels including sparsematrix multiple-vector operations and appropriate data structures for all relevant compute architectures provided the foundation for hardware-and energy-efficient large-scale computations.
Building on these high-performance building blocks, one focus of the algorithm layer was put on the block formulation of Jacobi-Davidson [64] and filter diagonalization [56] methods, the hardware efficiency of preconditioners [46][47][48], and the development of hardware-aware coloring schemes [1].In terms of scalability, the project has investigated new contour-based integration eigensolvers [23,24] that can exploit additional parallelism layers beyond the usual data parallelism.The solvers developed in ESSEX can tackle standard, generalized, and nonlinear eigenvalue problems and may also be used to extract large bulks of extremal and inner eigenvalues.
The applications layer applies the algorithms and building blocks to deliver scalable solutions for topical quantum materials like graphene, topological insulators, or superconductors, and nonlinear dynamical systems like reaction-diffusion systems.A key issue for large-scale simulations is the scalable (in terms of size and parallelism) generation of the sparse matrix representing the model Hamiltonian.Our matrix generation framework ScaMaC can be integrated into application code to allow the on-the-fly, in-place construction of the sparse matrix.Beyond the ESSEX application fields, matrices from many other relevant areas can be produced by ScaMaC.
The PHIST library [78] is the sustainable outcome of the performance-centric efforts in ESSEX.It is built on a rigorous software and performance engineering process, comprises diverse solver components, and supports multiple backends (e.g., Trilinos, PETSc, ESSEX kernels).It also interfaces to multiple languages such as C, C++, Fortran 2003, and Python.The CRAFT library [73] provides user-friendly access to fault tolerance via checkpoint/restart and automatic recovery for iterative codes using standard C++.
This review focuses on important developments in ESSEX-II.After presenting a brief overview of the most relevant achievements in the first project phase ESSEX-I in Sects. 2 and 3 details algorithmic developments in ESSEX-II, notably with respect to preconditioners and projection-based methods for obtaining inner eigenvalues.Moreover, we present the RACE (Recursive Algebraic Coloring Engine) method, which delivers hardware-efficient graph colorings for parallelization of algorithms and kernels with data dependencies.In Sect. 4 we showcase performance and parallel efficiency numbers for library components developed in ESSEX-II that are of paramount importance for the application work packages: GPGPU-based tall and skinny matrix-matrix multiplication and the computation of inner eigenvalues using polynomial filter techniques.Section 5 describes the software packages that were developed to a usable and sustainable state, together with their areas of applicability.In Sect.6 we show application results from the important areas of quantum physics and nonlinear dynamical systems.Finally, in Sect.7 we highlight the collaborations sparked and supported by SPPEXA through the ESSEX-II project.

Summary of the ESSEX-I Software Structure
The Exascale-enabled Sparse Solver Repository (ESSR) was developed along the requirements of the algorithms and applications under investigation in ESSEX.It was not intended as a full-fledged replacement of existing libraries like Trilinos 5 [33], but rather as a toolbox that can supply developers with blueprints as a starting point for their own developments.In ESSEX-I, the foundations for a sustainable software framework were laid.See Sect. 5 for developments in ESSEX-II.
The initial version of the ESSR [77] comprised four components: • GHOST (General, Hybrid and Optimized Sparse Toolkit) [55], a library of basic sparse and dense linear algebra building blocks that are not available in this form in other software packages.The development of GHOST was strictly guided by performance engineering techniques; implementations of standard kernels such as sparse matrix-vector multiplication (spMVM) and sparse matrix-multiple-vector multiplication (spMMVM) as well as tailor-made fused kernels, for instance those employed in the Kernel Polynomial Method (KPM) [81], were modeled using the roofline model.GHOST supports, by design, strongly heterogeneous environments using the MPI+X approach.See [51] for a comprehensive overview of GHOST and its building blocks.• ESSEX-Physics, a collection of prototype implementations of polynomial eigensolvers such as the KPM and Chebyshev Filter Diagonalization (ChebFD).These were implemented on top of GHOST using tailored kernels and were shown to perform well on heterogeneous CPU-GPU systems [53]. 5https://trilinos.org/.
• PHIST (Pipelined Hybrid-parallel Iterative Solver Toolkit), which comprises Jacobi-Davidson type eigensolvers and Krylov methods for linear systems.One important component is a test framework that allows for continuous integration (CI) throughout the development cycle.PHIST can not only use plain GHOST as its basic linear algebra layer; it is also equipped with fallback kernel implementations and adapters for the Trilinos and Anasazi libraries.A major achievement in the development of PHIST was an efficient block Jacobi-Davidson eigenvalue solver, which could be shown to have significant performance advantages over nonblocked versions when using optimized building blocks from GHOST [64].• BEAST (Beyond fEAST), which implements innovative projection-based eigensolvers motivated by the contour integration-based FEAST method [23].The ESSEX-I project has contributed to improving FEAST in two ways: by proposing techniques for solving or avoiding the linear systems that arise, and by improving robustness and performance of the algorithmic scheme.
A pivotal choice for any sparse algorithm implementation is the sparse matrix storage format.In order to avoid data conversion and the need to support multiple hardware-specific formats in a single code, we developed the SELL-Cσ format [52].It shows competitive performance on the dominating processor architectures in HPC: standard multicore server CPUs with short-vector single instruction multiple data (SIMD) capabilities, general-purpose graphics processing units (GPGPUs), and many-core designs with rather weak cores such as the Intel Xeon Phi.SELL-C-σ circumvents the performance penalties of matrices with few nonzero entries per row on architectures on which SIMD vectorization is a key element for performance even with memory-bound workloads.
In order to convert a sparse matrix to SELL-C-σ , its rows are first sorted according to their respective numbers of nonzeros.This sorting is performed across row blocks of length σ .After that, the matrix is cut into row blocks of length C. Within each block, rows are padded with zeros to equal length and then stored in column-major order.See Fig. 1 for visualizations of SELL-C-σ with C = 6 and σ ∈ {1, 12, 24}.Incidentally, known and popular formats can be recovered as corner cases: SELL-1-1 is the well-known CSR storage format and SELL-N-1 is ELLPACK.The particular choice of C and σ influences the performance of the spMVM operation; optimal values are typically matrix-and hardware-dependent.However, in practice one can usually find parameters that yield good performance across architectures for a particular matrix.A roofline performance model was constructed in [52] that sets an upper limit for the spMVM performance for any combination of matrix and architecture.This way, "bad" performance is easily identified.SELL-C-σ was quickly adopted by the community and is in use, in pure or adapted form, in many performance-oriented projects [5,6,28,60,80].

Algorithmic Developments
In this section we describe selected developments within ESSEX-II on the algorithmic level, in particular preconditioners for the solution of linear systems that occur in the eigensolvers, a versatile framework for computing inner eigenvalues, and a nonlinear eigensolver.We also cover a systematic comparison of contourbased methods.We close the section with the introduction of RACE, which is an algorithmic development for graph coloring guided by the constraints of hardware efficiency.

Preconditioners (ppOpen-SOL)
Two kinds of solvers have been developed: a preconditioner targeting the illconditioned large scale problems arising in the BEAST-C method (cf.Sect.3.2) and a multigrid solver targeting problems arising from finite difference discretizations of partial differential equations (PDEs).

Regularization
The BEAST-C method leads to a large number of ill-conditioned linear systems with complex diagonal shifts [24].Furthermore, in many of our quantum physics applications, the system matrices have small (and sometimes random) diagonal elements.In order to apply a classic incomplete Cholesky (IC) factorization preconditioner, we used two types of regularization to achieve robustness: a blocking technique (BIC) and an additional diagonal shift [47].Using this approach, we solved a set of 120 prototypical linear systems from this context (e.g., BEAST-C applied to quantum physics applications).Due to the complex shift, the system matrix is symmetric but not Hermitian.Hence we use an adaptation of the Conjugate Gradient (CG) method for complex symmetric matrices called COCG (conjugate orthogonal conjugate gradient [79]).
The blocking technique is a well-known approach for improving the convergence rate.In this study, we apply the technique not only for better convergence but also for more robustness.The diagonal entries in the target equations are small.By applying the blocking technique, the diagonal blocks to be inverted include larger off-diagonal entries.
The diagonal shifting is a direct measure for transforming the ill-conditioned matrices to be more diagonally dominant before performing the incomplete factorization.On the other hand, this may deteriorate the convergence of the overall method.We therefore investigate the best value for the diagonal shifting for our applications.
Figure 2 shows the effect of the regularized IC preconditioner with the COCG method.By using the diagonal shifted block IC-COCG (BIC-COCG), we solve all target linear systems.
Fig. 2 Effect of the regularized IC preconditioner with the COCG method.By using the diagonal shifted block IC-COCG (BIC-COCG), we can solve all test problems from our benchmark set

Hierarchical Parallel Reordering
In this section, we present scalability results for the BIC preconditioner parallelized by a hierarchical parallel graph coloring algorithm.This approach yields an almost constant convergence rate with respect to the number of compute nodes, and good parallel performance.
Node-wise multi-coloring (with domain decomposition between nodes) is widely used for parallelizing IC preconditioners on clusters of shared memory CPUs.Such "localized" multi-coloring leads to a loss of robustness of the regularized IC-COCG method, and the convergence rate decreases at high levels of parallelism.To solve this problem, we parallelize the block IC preconditioner for the hybrid-parallel cluster system.In addition, we proposed the hierarchical parallelization for the multi-coloring algorithms [46].This versatile scheme allows us to parallelize almost any multi-coloring algorithm.
Figure 3 shows the number of iterations and computational time of the BIC-COCG method on the Oakleaf-FX cluster, using up to 4,800 nodes.The benchmark matrix is the Hamiltonian of a graphene sheet simulation with more than 500 million linear equations, for which interior eigenvalues are of interest [24].Hierarchical parallelization yields almost constant convergence with respect to the number of nodes.The computational time with 4,600 nodes is 30 times smaller than with 128 nodes, amounting to a parallel efficiency of 83.5% if the 128-node case is taken as the baseline.

Multiplicative Schwarz-Type Block Red-Black Gauß-Seidel Smoother
Multigrid methods are among the most useful preconditioners for elliptic PDEs.
In [45] we proposed a multiplicative Schwarz block red/black Gauß-Seidel (MS- The unknowns are divided into blocks so that the amount of data for processing each block fits into the cache, and α Gauß-Seidel iterations are applied to the block per smoother step.The computational cost for the additional iterations is much lower than for the first iteration because of data locality. Figure 4 shows the effect of the MS-BRB-GS(α) smoother on a single node of the ITO system (Intel Xeon Gold 6154 (Skylake-SP) Cluster at Kyushu University).By increasing the number of both pre-and post-smoothing steps, the number of iterations is decreased.In the best case, MS-BRB-GS is 1.64× faster than BRB-GS.

The BEAST Framework for Interior Definite Generalized Eigenproblems
The BEAST framework targets the solution of interior definite eigenproblems In the following we highlight some of BEAST's algorithmic features, skipping other topics such as locking converged eigenpairs, adjusting the dimension of the subspace, and others.

Projector Types
BEAST provides three variants of approximate projectors.First, polynomial approximation (BEAST-P) using Chebyshev polynomials, which only requires matrix vector multiplications but is restricted to standard eigenproblems.Second, Cauchy integral-based contour integration (BEAST-C), as in the FEAST method [63].As a third method, an iterative implementation of the Sakurai-Sugiura method [65] is available (BEAST-M), which shares algorithmic similarities with FEAST.In the following we briefly elaborate on the algorithmic ideas.
• In BEAST-P, we have U = p(A) • Y with a polynomial p(z) = d k=0 c k T k (z) of suitable degree d.Here, T k denotes the kth Chebyshev polynomial, Due to the use of the T k , this method is also known as Chebyshev filter diagonalization.In addition to well-known methods for computing the coefficients c k [18,62], BEAST also provides the option of using new, improved coefficients [25].Their computation depends on two parameters, μ and σ , and for suitable combinations of these, the filtering quality of the polynomial can be improved significantly; see Fig. 5, which shows the "gain," i.e., the reduction of the width of those λ values outside the search interval, for which a damping of corresponding eigenvectors by at least a factor 100 cannot be guaranteed.For some combinations (σ, μ), marked red in the picture, this "no guarantee" area can be reduced by a factor of more than 2, which in turn allows using lower-degree polynomials to achieve comparable overall convergence.A parallelized method for finding suitable parameter combinations and computing the c k is included with BEAST.
• In BEAST-C, the exact projection (integration is over a contour in the complex plane that encloses the eigenvalues λ ∈ [λ, λ], but no others) is approximated using an N-point quadrature rule, leading to N linear systems, where the number of right-hand sides (RHS) corresponds to the dimension of the current subspace U (and Y ).• BEAST-M is also based on contour integration, but moments are used to reduce the number of RHS in the linear systems.Taking M moments, we have and thus an M times smaller number of RHS (dimension of Y ) is sufficient to achieve the same dimension of U .
The linear systems in the contour-based schemes may be ill-conditioned if the integration points z j are close to the spectrum (this happens, e.g., for narrow search intervals [λ, λ]); cf. also Sects.3.1 and 3.4 for approaches to address this issue.

Flexibility, Adaptivity and Auto-Tuning
The BEAST framework provides flexibility at the algorithmic, parameter, and working precision levels, which we describe in detail in the following.

Algorithmic Level
The projector can be chosen from the three types described above, and the type may even be changed between iterations.In particular, an innovative subspaceiterative version of Sakurai-Sugiura methods (SSM) has been investigated for possible cost savings in the solution of linear systems via a limited subspace size and the overall reduction of number of right hand sides over iterations by using moments.Given, however, the potentially reduced convergence threshold with a constrained subspace size, we support switching from the multi-moment method, BEAST-M, to a single-moment method, BEAST-C.The efficiency, robustness, and accuracy of this approach in comparison with traditional SSM and FEAST has been explored [35].We further studied this scheme along with another performance-based implementation of SSM, z-PARES [65,67].These investigations considered the scaling and computational cost of the libraries as well as heuristics for parameter choice, in particular with respect to the number of quadrature nodes.We observed that the scaling behavior improved when the number of quadrature nodes increased, as seen in Fig. 6.As the linear systems solved at each quadrature node are independent and the quality of numerical integration improves with increased quadrature degree, exploiting this property makes sense, particularly within the context of exascale computations.However, it is a slightly surprising result, as previous experiments with FEAST showed diminishing returns for convergence with increased quadrature degree [23], something we do not observe here.

Parameter Level
In addition to the projector type, several algorithmic parameters determine the efficiency of the overall method, most notably the dimension of the subspace and the degree of the polynomial (BEAST-P) or the number of integration nodes (BEAST-C and BEAST-M).
With certain assumptions on the overall distribution of the eigenvalues, clear recommendations for optimum subspace size (as a multiple of the number of expected eigenvalues) and the degree can be given, in the sense that overall work N refers to the number of quadrature nodes along a circular contour, NP to the number of processes is minimized.For more details, together with a description of a performance-tuned kernel for the evaluation of p(A) • Y , the reader is referred to [62].
If such information is not available, or for the contour integration-type projectors, a heuristic has been developed that automatically adjusts the degree (or number of integration nodes) during successive iterations in order to achieve a damping of the unwanted components by a factor of 100 per iteration, which leads to close-tooptimum overall effort; cf.[26].

Working Precision Level
Given the iterative nature of BEAST, with one iteration being comparatively expensive, the possibility to reduce the cost of at least some of these iterations is attractive.We have observed that before a given residual tolerance is surpassed, systematic errors in the computation of the projector and other operations do not impair convergence speed per se, but impose a limit on what residual can be reached before progress stagnates.One such systematic error is the finite accuracy of floating-point computations, which typically are available in single and double precision.In the light of the aforementioned behavior, it seems natural to perform initial iterations in single precision and thereby save on computation time before a switch to double precision becomes inevitable; cf.Fig. 7. Therefore, mixed precision has been implemented in all BEAST schemes mentioned above, allowing an adaptive strategy to automatically switch from single to double precision after a given residual tolerance is reached.A comprehensive description and results are presented in [4].These results and our initial investigations also suggest that increased precision beyond double precision (i.e., quad precision) will have no benefit for the convergence rate until a certain double precision specific threshold is reached; convergence beyond this point would require all operations to be carried out with increased precision.

Levels of Parallelism
The BEAST framework exploits multiple levels of parallelism using an MPI+X paradigm.We rely on the GHOST and PHIST libraries for efficient sparse matrix/dense vector storage and computation; cf.Sect. 5.The operations implemented therein are themselves hybrid parallel and constitute the lowest level of parallelism in BEAST.Additional levels are addressed by parallelizing over blocks of vectors in Y and, for BEAST-C and BEAST-M, over integration nodes during the application of the approximate projector.A final level is added by exploiting the ability of the method to subdivide the search interval [λ, λ] and to process the subintervals independently and in parallel.Making use of these properties, however, may lead to non-orthogonal eigenvectors, which necessitates postprocessing as explained in the following.

A Posteriori Cross-Interval Orthogonalization
Rayleigh-Ritz-based subspace iteration algorithms naturally produce a Borthogonal set of eigenvectors X, i.e., orth(X) is small, where orth(X) = max orth(x i , x j )|i = j with orth(x, y) = y, x x y .
By contrast, the orthogonality orth(X, Y ) = max orth(x i , y j ) between two or more independently computed sets of eigenvectors may suffer if the distance between the involved eigenvalues is small [49,50].Simultaneous reorthogonalization of evolving approximate eigenvectors during subspace iteration has proven ineffective unless the vectors have advanced reasonably far.A large scale re-orthogonalization of finished eigenvector blocks, on the other hand, requires a careful choice of methodology in order to not diminish the quality of the previously established residual.
Orthogonalization of multiple vector blocks implies Gram-Schmidt style propagation of orthogonality, assuming orth(X, Y ) can be arbitrarily poor.In practice, the independently computed eigenvectors will exhibit multiple grades of orthogonality, but rarely will there be no orthogonality (in the sense above) at all.This, in turn, allows for the use of less strict orthogonalization methods.While, in theory, the orthogonalization of p blocks requires at least p(p − 1)/2 + (p − 1) block-block or intra-block orthogonalizations and ensures global orthogonality, an iterative scheme allows for more educated choices on the ordering of orthogonalizations in order to reduce losses in residual and improve the communication pattern, eliminating the need for broadcasts of vector blocks at the cost of additional orthogonalization operations in the form of multiple sweeps.In practice, very few sweeps (∼2) are sufficient in most cases.
Every block-block orthogonalization X = X − Y (Y H BX) disturbs the orthogonality orth(X) of the modified block, as well as its residual.Local re-orthogonalization of X disturbs the residual further.We have identified orthogonalization patterns and selected orthogonalization algorithms that reduce the loss of residual accuracy to a degree that essentially eliminates the need for additional post-iteration.
The implementation of an all-to-all interaction of many participating vector blocks can be performed in multiple ways with different requirements regarding storage, communication, runtime, and with different implications on accuracy and loss of residual.Among several such strategies and algorithms that have been implemented and tested, the most promising is a purely iterative scheme, both for global and local orthogonalization operations.It is based on a comparison of interval properties, most notably the achieved residual from the subspace iteration.We are continuing to explore the possibility to detect certain orthogonalizations as unnecessary without computing the associated inner products in order to further reduce the workload without sacrificing orthogonality.

Robustness and Resilience
In the advent of large scale HPC clusters, hardware faults, both detectable an undetectable, have to be expected.
Detectable hardware faults, e.g., the outage of a component that violently halts execution, can typically only be mitigated by frequent on-the-fly storage of the most vital information.In the case of subspace iteration, as is used in BEAST, almost all required information for being able to resume computation is encoded in the iterated subspace basis in form of the approximate eigenvectors, besides runtime information about the general program flow.Relying on the CRAFT library [73], a per-iteration checkpointing mechanism has been implemented in BEAST.
Additionally, for also being able to react to "silent" computation errors that merely distort the results but do not halt execution, the most expensive operation (application of the approximate projector) has been augmented to monitor the sanity of the results.This can be done in two ways: A checksum-style entrainment of additional vectors, linear combinations of the right-hand sides, can be checked during and after the application of the projector to detect errors and allow for the re-computation of the incorrect parts.The comparison of approximate filter values obtained from the computed basis and the expected values obtained from the scalar representation of the filter function, on the other hand, gives an additional a posteriori test for the overall plausibility of the basis.
Practical tests have shown that small distortions of the subspace basis have not enough impact on the overall process in order to justify expensive measures.If the error is not recurring, just continuing the subspace iteration is often the best and most cost-efficient option.This is particularly true in early iterations, where small errors have no effect at all.

Relationship Among Contour Integral-Based Eigensolvers
The complex moment-based eigensolvers such as the Sakurai-Sugiura method can be regarded as projection methods using a subspace constructed by the contour integral The property of the subspace is well analyzed by using a filter function which approximates a band-pass filter for the target region where the wanted eigenvalues are located.Using the filter function, error analyses of the complex moment-based eigensolvers were shown in [30,41,42,66,76].By using the results of the error analyses, an error resilience technique and an accuracy deterioration technique have also been given in [32,43].
The relationship between typical complex moment-based eigensolvers was also analyzed in [42] focusing on the subspace.The block SS-RR method [36] and the FEAST algorithm [76] are projection methods for solving the target generalized eigenvalue problem, whereas the block SS-Hankel method [37], Beyn [12], the block SS-Arnoldi methods [40] and its improvements [38] are projection methods for solving an implicitly constructed standard eigenvalue problem; see [42] for details. Figure 8 shows a map of the relationships among the contour integral-based eigensolvers.

Extension to Nonlinear Eigenvalue Problems
The complex moment-based eigensolvers were extended to nonlinear eigenvalue problems (NEPs): where the matrix-valued function T : → C n×n is holomorphic in an open domain .The projection for a nonlinear matrix function T (λ) is given by This projection is approximated by The block SS-Hankel [7,8], block SS-RR [83], and block SS-CAA methods [39] are simple extensions of the GEP solvers.A technique for improving the numerical stability of the block SS-RR method for NEP was developed in [15,16].
Beyn proposed a method using Keldysh's theorem and the singular value decomposition [12].Van Barel and Kravanja proposed an improvement of the Beyn method using the canonical polyadic (CP) decomposition [10].

Recursive Algebraic Coloring Engine (RACE)
The standard approach to solve the ill-conditioned linear systems arising in BEAST-C or FEAST is to use direct solvers.However, in [24] it was shown that the Kaczmarz iterative solver accelerated by a Conjugate Gradient (CG) method (the so-called CGMN solver [29]) is a robust alternative to direct solvers.Standard multicoloring (MC) was used in [29] for the parallelization of the CGMN kernels.After analyzing the shortcomings of this strategy in view of hardware efficiency, we developed in collaboration with the EXASTEEL-II project the Recursive Algebraic Coloring Engine (RACE) [1].It is an alternative to the well-known MC and algebraic block multicoloring (ABMC) algorithms [44], which have the problem that their matrix reordering can adversely affect data access locality.RACE aims at improving data locality, reducing synchronization, and generating sufficient parallelism while still retaining simple matrix storage formats such as compressed row storage (CRS).We further identified distance-2 coloring of the underlying graph as an opportunity for parallelization of the symmetric spMVM (SymmSpMV) kernel.
RACE is a sequential, recursive, level-based algorithm that is applicable to general distance-k dependencies.It is currently limited to matrices with symmetric structure (undirected graph), but possibly nonsymmetric entries.The algorithm comprises four steps: level construction, permutation, distance-k coloring, and load balancing.If these steps do not generate sufficient parallelism, recursion on subgraphs can be applied.Using RACE implies a pre-processing and a processing phase.In pre-processing, the user supplies the matrix, the kernel requirements (e.g., distance-1 or distance-2) and hardware settings (number of threads, affinity strategy).The library generates a permutation and stores the recursive coloring information in a level tree.It also creates a pool of pinned threads to be used later.In the processing phase, the user provides a sequential kernel function which the library executes in parallel as a callback using the thread pool.
Figure 9 shows the performance of SymmSpMV on a 24-core Intel Xeon Skylake CPU for a range of sparse symmetric matrices.In Fig. 9a we compare RACE against Intel's implementation in the MKL library, and with roofline limits obtained via bandwidth measurements using array copy and read-only kernels, respectively.RACE outperforms MKL by far.In Fig. 9b we compare against standard multicoloring (MC) and algebraic block multicoloring (ABMC).The advantage of RACE is especially pronounced with large matrices, where data traffic and locality of access is pivotal.One has to be aware that some algorithms may exhibit a change in convergence behavior due to the reordering.This has to be taken into account when benchmarking whole program performance instead of kernels.Details can be found in [1].
In order to show the advantages of RACE in the context of a relevant algorithm, we chose FEAST [63] for computing inner eigenvalues.The hot spot of the algorithm (more than 95%) is a solver for shifted linear systems (A−σ I = b).These systems are, however, highly ill-conditioned, posing severe convergence problems for most linear iterative solvers.We use the FEAST implementation of Intel MKL, which by default employs the PARDISO direct solver [69], but its Reverse Communication Interface (RCI) allows us to plug our CGMN implementation instead.In the following experiment we find ten inner eigenvalues of a simple discrete Laplacian matrix to an accuracy of 10 −8 .Figure 10 shows the measured time and memory footprint of the default MKL version (using PARDISO) and the CGMN versions parallelized using both RACE and ABMC for different matrix sizes.ABMC is a factor of 4× slower than RACE.The time required by the default MKL with PARDISO is smaller than with CGMN using RACE for small sizes; however, the gap gets smaller as the size grows due to the direct solvers having a higher time complexity (here ≈ O(n 2 )) compared to iterative methods (≈ O(n 1.56 )).Moreover, the direct solver requires more memory, and the memory requirement grows much faster (see Fig. 10b) than with CGMN.In our experiment the direct solver ran out of memory at problem sizes beyond 140 3 , while CGMN using RACE used less than 10% of space at this point.Thus, CGMN with RACE can solve much larger problems compared to direct solvers, which is a major advantage in fields like quantum physics.

Hardware Efficiency and Scalability
In this section we showcase performance and parallel efficiency numbers for library components developed in ESSEX-II that are of paramount importance for the application work packages: GPGPU-based tall and skinny matrix-matrix multiplication and the computation of inner eigenvalues using polynomial filter techniques.

Tall and Skinny Matrix-Matrix Multiplication (TSMM) on GPGPUs
Orthogonalization algorithms frequently require the multiplication of matrices that are strongly nonsquare.Vendor-supplied optimized BLAS libraries often yield suboptimal performance in this case."Sub-optimal" is a well-defined term here since the multiplication of an M×K matrix A with a K×N matrix B with K M, N and small M, N is a memory-bound operation: At M = N, its computational intensity is just M/8 flop/byte.In ESSEX-I, efficient implementations of TSMM on multicore CPUs were developed [51].
The naive roofline model predicts memory-bound execution for M 64 on a modern Volta-class GPGPU.See Fig. 11 for a comparison of optimal (roofline) performance and measured performance for TSMM on an Nvidia Tesla V100 GPGPU using the cuBLAS library. 6We have developed an implementation of TSMM for GPGPUs [20], investigating various optimization techniques such as different thread mappings, overlapping long-latency loads with computation via leapfrogging 7 and unrolling, options for global reductions, and register tiling.Due Fig. 12 Best achieved performance for each matrix size with M = N in comparison with the roofline limit, cuBLAS and CUTLASS, with K = 2 23 .(From [20]) to the large and multi-dimensional parameter space, the kernel code is generated using a python script.
Figure 12 shows a comparison between our best implementations obtained via parameter search (labeled "leap frog" and "no leap frog," respectively) with cuBLAS and CUTLASS, 8 which is a collection of CUDA C++ template abstractions for high-performance matrix multiplications.Up to M = N = 36, our implementation stays within 95% of the bandwidth limit.Although the performance 8 https://github.com/NVIDIA/cutlass(May 2019).levels off at larger M, N, which is due to insufficient memory parallelism, it is still significantly better than with cuBLAS or CUTLASS.

Node-Level Performance
Single-device benchmark tests for BEAST-P were performed on an Intel Knights Landing (KNL), an Nvidia Tesla P100, and an Nvidia Tesla V100 accelerator, comparing implementations based on vendor libraries (MKL and cuBLAS/cuSPARSE, respectively) with two versions based on GHOST: one with and one without tailored fused kernels.The GPGPUs showed performance levels expected from a bandwidthlimited code, while on KNL the bottleneck was located in the core (see Fig. 13a).Overall, the concept of fused optimized kernels provided speedups of up to 2× compared to baseline versions.Details can be found in [56].

Massively Parallel Performance
Scaling tests for BEAST-P were performed on the "Oakforest-PACS" (OFP) at the University of Tokyo, "Piz Daint" at CSCS in Lugano, and on the "SuperMUC-NG" (SNG) at Leibniz Supercomputing Centre (LRZ) in Garching. 9While the OFP nodes comprise Intel "Knights Landing" (KNL) many-core CPUs, SNG has CPU-only dual-socket nodes with Intel Skylake-SP, and Piz Daint is equipped with single-socket Xeon "Haswell" nodes, each of which has an Nvidia Tesla P100 accelerator attached.Weak and strong scaling tests were done with topological insulator (TI) matrices generated by the ScaMaC library.Flops were calculated for the computation of the approximate eigenspace, U , averaged over the four BEAST iterations it took to find the 148 eigenvalues in each interval to a tolerance of 1 × 10 −10 .The subspace contained 256 columns, and spMMVs were performed in blocks of size 32 for best performance.Optimized coefficients [25] were used for the Chebyshev polynomial approximation, resulting in a lower overall required polynomial degree.Weak and strong scaling results are shown in Fig. 14a through d.
OFP and SNG show similar weak scaling efficiency due to comparable singlenode performance and network characteristics.Piz Daint, owing to its superior single-node performance of beyond 400 Gflop/s, achieves only 60% of parallel efficiency at 2048 nodes.A peculiar observation was made on the CPU-only SNG system: Although the code runs fastest with pure OpenMP on a single node (223 Gflop/s), scaled performance was observed to be better with one MPI process per socket.The ideal scaling and efficiency numbers in Fig. 14a-c use the best value on the smallest number of nodes in the set as a reference.The largest matrix on SNG had 6.6 × 10 9 rows.

Scalable and Sustainable Software
It was a central goal of the ESSEX-II project to consolidate our software efforts and provide a library of solvers for sparse eigenvalue problems on extreme-scale HPC systems.This section gives an overview of the status of our software, most of which is now publicly available under a three-clause BSD license.Many of the efforts have been integrated in the PHIST library so that they can easily be used together, and we made part of the software available in larger contexts like Spack [27] and the extreme-scale scientific software development kit xSDK [11].The xSDK is an effort to define common standards for high-performance, scientific software in terms of software engineering and interoperability.Dashed lines denote ideal scaling with respect to the smallest number of nodes in the set.(a) Weak scaling of BEAST-P on OFP for problems of size 2 20 (2 nodes) to 2 30 (2048 nodes, about one quarter of the full machine).(b) Weak scaling of BEAST-P on Piz Daint for problems of size 2 21  (1 node) to 2 32 (2048 nodes).(From [56]).(c) Weak scaling of BEAST-P on SNG for problems of size 2 21 (1 node) to 1.53 × 2 32 (3136 nodes, about half of the full machine).(d) Strong scaling of BEAST-P on SNG for problems of size 2 28 (crosses) and 2 30 (triangles) The current status of the software developed in the ESSEX-II project is summarized as follows.
• BEAST is available via bitbucket, 10 and can be compiled either using the PHIST kernel interface or the GHOST library directly.The former allows using it with any backend supported by PHIST.• CRAFT is available stand-alone 11 or (in a fixed version) as part of PHIST.
• ScaMaC is available stand-alone 12 or (in a fixed version) as part of PHIST.
• GHOST is available via bitbucket. 13The functionality which is required to provide the PHIST interface can be tested via PHIST.Achieving full (or even substantial) test coverage of the GHOST-functionality would require a very large number of tests (in addition to what the PHIST interface provides, GHOST allows mixing data types, and it uses automatic code generation, which leads to an exponentially growing number of possible code paths with every new kernel, supported processor and data type).It is, however, possible to create a basic GHOST installation via the Spack package manager (since March 2018, commit bcde376).• PHIST is available via bitbucket 14 and Spack (since commit 2e4378b).
Furthermore, PHIST 1.7.5 is part of xSDK 0.4.0.The version distributed with the xSDK is restricted to use the Tpetra kernels to maximize the interoperability of the package.

PHIST and the Block-ILU
In ESSEX-I we addressed mostly node-level performance [77] on multi-core CPUs.
The main publication of ESSEX-II concerning the PHIST library [78] presents performance results for the block Jacobi-Davidson QR (BJDQR) solver on various platforms, including recent CPUs, many-core processors and GPUs.It was also shown in this work that the block variant has a clear performance advantage over the single-vector algorithm in the strong scaling limit.The reason is that, while the number of matrix-vector multiplications increases with the block size (see also [64]), the total number of reductions decreases.In order to demonstrate the performance portability of PHIST, we show in Fig. 15 a weak scaling experiment on the recent SuperMUC-NG machine.
For the block size 4, we roughly match the performance it achieves in the memory-bounded HPCCG benchmark (207 TFlop/s), 15 but using only half of the machine.This gives a clear indication that our node-level performance engineering and multi-node implementation are highly successful: after all, we do not optimize for the specific operator application (a simple structured grid, 3D Laplace operator), which the HPCCG code does.On the other hand, we have an increased computational intensity for some of the operations due to the blocking, which increases the performance over a single-vector CG solver.The single-vector BJDQR solver achieves 98 TFlop/s on half of the machine.

Integration of the Block-ILU Preconditioning Technique
Initial steps have been taken to make the Block-ILU preconditioner (cf.Sect.3.1) available via the PHIST preconditioning interface.At the time of writing, there is an experimental implementation of a block CRS sparse matrix format in the PHIST builtin kernel library, including parallel conversion and matrix-vector product routines and the possibility to construct and apply the block Cholesky preconditioner.Furthermore, the interfaces necessary to allow using the preconditioner within the BJDQR eigensolver have been implemented.These features are available for experimenting in a branch of the PHIST git repository because they do not yet meet the high demands on maintainability (especially unit testing) and documentation of a publicly available library.Integration of the method with the BEAST eigensolver is not yet possible because the builtin kernel library does not support complex arithmetic.As mentioned in Sect.3.2, the complex version will be integrated directly into the BEAST software, instead.

BEAST
BEAST combines implementations of spectral filtering methods for Rayleigh-Ritz type subspace iteration in a generalized framework to provide facilities for improving performance and robustness.The algorithmic foundation allows for the solution of interior Hermitian definite eigenproblems of standard and generalized form via an iterative eigensolver, unveiling all eigenpairs in one or many specified intervals.The software is designed as hybrid parallel library, written in C/C++, and relying on GHOST and PHIST to provide basic operations, parallelism, and data types.Beyond the excellent scalability of the underlying kernel libraries, multiple additional levels of parallelism allow for computing larger portions of the spectrum and/or utilizing a larger number of computing cores.The inherent ability of the underlying algorithm to compute separate intervals independently offers wide potential but requires careful handling of cross-interval interactions to ensure the desired quality of results, which is well supported by BEAST.
The BEAST library interface comes in variations for the common floating point formats (real and complex, single and double precision) for standard and generalized eigenproblems.Additionally, the software offers the possibility to switch precisions on-the-fly, from single to double precision, in order to further improve performance.While BEAST offers an algorithm for standard eigenproblems that completely bypasses the need for linear system solves, other setups typically require a suitable linear solver.Besides a builtin parallel sparse direct solver for banded systems, BEAST includes interfaces to MUMPS and Strumpack, as well as a flexible callback-driven interface for the inclusion of arbitrary linear solvers.It also interfaces with CRAFT and ScaMaC, which provide fault tolerance and dynamic matrix generation, respectively.While working out of the box for many problems, BEAST offers a vast amount of options to tweak the software for the specific problem at hand.A builtin command line parser allows for easy modification.The included application bundles the several capabilities of BEAST in form of a standalone tool that reads or generates matrices and solves the specified eigenproblem.As such, it acts as comprehensive example for the usage of BEAST.
The library is still in a development state, and interface and option sets may change.A more comprehensive overview over a selection of features is provided in Sect.3.2.

CRAFT
The CRAFT library [73] covers two essential aspects of fault tolerance namely communication, and data recovery of an MPI application in case of process-failures.
In the Checkpoint/Restart part of the library, it provides an easier and extensible interface for making application-level checkpoint/restart.A CRAFT checkpoint can be defined simply by defining a Checkpoint object and adding the restartrelevant data in it, as shown in Listing 1.By default, the Checkpoint::add() function supports the most frequently used data formats, e.g., "plain old data" (POD), i.e., int, double, float, etc., POD 1D-and 2D-arrays, MPI data-types, etc..However, it can be easily extended to support any user defined data-types.The Checkpoint::read(), write() and update() methods can then be used to read/write all added checkpoint's data.The library supports asynchronous-1 # i n c l u d e <mpi .h> 2 # i n c l u d e < c r a f t .h> 3 i n t main ( i n t argc , c h a r * argv [ ] ) { 4 . . .Listing 1 A toy-code that demonstrates the simplicity of CRAFT's checkpoint/restart and automatic fault tolerace features in a typical iterative-style scientific application checkpointing as well as node-level checkpointing using the SCR library [68].Moreover it supports multi-staged, nested-, and signal-checkpointing.
The Automatic Fault Tolerance (AFT) part of CRAFT provides an easier interface for a dynamic process-failure recovery and management.CRAFT uses the ULFM-MPI implementation for process-failure detection, propagation, and communication recovery procedures, however it considerably reduces the user's effort by hiding these details behind AFT_BEGIN() and AFT_END() functions as shown in Listing 1.After a process failure, the library recovers the broken communicator (shrinking or non-shrinking by process-spawning), and returns the control back to the program at AFT_BEGIN(), where the data can be recovered.Both of these CRAFT functionalities are designed to complement each other, however they can be used independently as well.For detailed explanation of the features included in CRAFT, check [73].Moreover, the library is available at [72].

CRAFT Benchmark Application
Within the scope of ESSEX, we have integrated CRAFT in the GHOST and PHIST libraries, and the BEAST algorithm.
Figure 16 shows a benchmark comparing the overhead of three different checkpointing strategies for the Lanczos algorithm (GHOST-based eigensolver), Jacobi-Davidson (PHIST-based eigensolver), and the BEAST algorithm.The important Fig. 16 CRAFT checkpointing overhead comparison for the Lanczos, Jacobi-Davidson, and BEAST eigenvalue solvers using three checkpointing methods of CRAFT, namely, node-level checkpointing with SCR, asynchronous PFS, and synchronous PFS checkpoints.The overhead for each checkpoint case is shown as a percentage.(Number of nodes=128, number of processes=256, Intel MPI) parameters for these benchmarks are listed in Table 1.The benchmark shows that the node-level and asynchronous checkpointing significantly reduces the checkpoint overhead despite a very high checkpoint frequency.
The benchmark presented in Fig. 17 demonstrates the overhead caused by checkpoint/restart as well as by the communication recovery after process failures for the Lanczos application.The first two bars, namely 'No CP Intel MPI' and 'No CP ULFM-MPI' show the runtime between non-fault-tolerant (Intel-MPI) vs. a fault-tolerant MPI implementation (ULFM-MPI), and creates a baseline for ULFM-MPI implementation without any failures.The next two groups of bars show the application runtime with 0-,1-, and 2-failures with checkpoints taken on PFS-and node-level.The failures are triggered at the mid-point of two successive checkpoints from within the application to have a deterministic re-computation time, where each failure simulates a complete node-crash (2 simultaneous process failures) and recovery is performed in a non-shrinking fashion on spare nodes.The largest contribution to the overhead is caused by the re-computation part, whereas the communication repair overhead takes an average of ≈ 2.6 s only.
Besides ESSEX, CRAFT has been utilized in [22] to create a process-level fault tolerant FEM code based on the shrinking recovery style.Moreover, CRAFT has been recently integrated in the EXASTEEL [21] project.

ScaMaC
Sparse matrices are central objects in the ESSEX project because of its focus on large-scale numerical linear algebra problems.A sparse matrix, whether derived from the Hamiltonian of a quantum mechanical system, from the Laplacian in a partial differential equation, or simply given as an abstract entity with unknown properties, defines a problem to be solved.The solution may then consist of a set of eigenvalues and eigenvectors computed with the BEAST or Jacobi-Davidson algorithms or, more moderately, of an estimate of some matrix norm or the spectral radius.
Testing and benchmarking of linear algebra algorithms, but also of computational kernels such as spMVM, requires matrices of different type and different size.Standard collections such as the Matrix Market [58] or Florida Sparse Matrix Collection [17] cover a wide range of examples, but mainly provide matrices of fixed moderate size.As algorithms and implementations improve, such matrices become readily too small and limited to serve as realistic test and benchmark cases.
We therefore decided in the ESSEX project to establish a collection of scalable matrices-the ScaMaC.Every matrix in ScaMaC is parameterized by individual parameters that allow the user to scale up the matrix dimension and to modify other, for example spectral, properties of the matrix.ScaMaC includes simple test and benchmark matrices but also 'real-world' matrices from research studies and applications.A major goal of ScaMaC is to provide a flexible yet generic interface for matrix generation, together with the necessary infrastructure to allow for immediate access to the collection irrespective of the concrete usage case.
The ScaMaC approach to matrix generation is straightforward and simple: Matrices are generated row-by-row (or column-by-column).The entire complexity of the actual generation technique, which depends on the specific matrix example, is encapsulated in a ScamacGenerator type and hidden from the user.ScaMaC provides routines to create and destroy such a matrix generator, to query matrix parameters prior to the actual matrix generation, and to obtain each row of the matrix.The ScaMaC interface is entirely generic and identical for all matrices in the collection.
A minimal code example is given in Fig. 18.In this example, the matrix and its parameters are set by parsing an argument string of the form "MatrixName, parameter=...,..." in line 3, before all rows are generated in the loop in lines 12-17.As this examples shows, parallelization of matrix generation is not part of the ScaMaC, but lies within the responsibility of the calling program.All ScaMaC routines are thread-safe and can be embedded directly into MPI processes and OpenMP threads.This approach guarantees full flexibility for the user and is easily integrated into existing parallel matrix frameworks such as PETSc or Trilinos.Both BEAST and PHIST provide direct access to the ScaMaC, therefore freeing the user from any additional considerations when using ESSEX software.The key feature of ScaMaC is scalability, since the matrix rows (or columns) can be generated independently and in arbitrary order.For example at Oakforest-PACS (see Sect. 4.2.2), a Hubbard matrix (see below) with dimension ≥9 × 10 9 and ≥1.5 × 10 11 non-zeros is generated in less than a minute, using 2 10 MPI processes each of which generates an average of 1.5 × 10 5 rows per second.As explained, the task of efficiently storing or using the matrix is left to the calling program.
ScaMaC is accompanied by a small toolkit for exploration of the collection.The toolkit addresses some basic tasks such as querying matrix information or plotting the sparsity pattern, but is not intended to compete with production-level code or full-fledged solver libraries, which the ESSEX project provides with BEAST, GHOST, and PHIST.
At the moment, 16 the matrix generators included in ScaMaC strongly reflect our personal research interests in quantum physics, but the ScaMaC framework is entirely flexible and allows for easy inclusion of new types of matrices, provided that they can be generated in a scalable way.The next update (scheduled for 2020) will extend ScaMaC primarily with matrix generators for standard partial differ-ential equations, including stencil matrices and finite element discretizations for advection-diffusion and elasticity problems, wave propagation, and the Schrödinger equation.Additional examples that are well suited for scalable generation are regular and irregular graph Laplacians, which have gained renewed interest in the context of machine learning [14,59].
To obtain an idea of the 'real-world' application matrices already contained in ScaMaC, consider two examples: The celebrated Hubbard model of condensed matter physics (Hubbard) [34] and a theoretical model for excitons in the cuprous oxide from our own research in this field (Exciton) [2].These matrices appear as Hamiltonians in the Schrödinger equation, and thus are either symmetric real (Hubbard) or Hermitian complex (Exciton).The respective application requires a moderate number (typically, 10-1000) of extremal or interior eigenpairs, which is less than 0.1% of the spectrum.Other ScaMaC generators provide general (non-symmetric or non-Hermitian) matrices, with a variety of sparsity patterns, spectral properties, etc.All generators depend on a number of application-specific parameters, 17 which are partly listed in Table 2 for the Hubbard and Exciton generator.
For the Hubbard example, two parameters determine the matrix dimension and sparsity pattern: n_fermions gives the number electrons per spin orientation (up or down), n_sites the number of orbitals occupied by the electrons.In terms of these parameters, the matrix dimension is D = n_sites n_fermions 2 .This dependency results in the rapid growth of D shown in Table 3.In the physically very interesting case of half-filling (n_fermions = n_sites/2 = n) we have asymptotically D 2 n / √ (π/2)n, that is, exponential growth of D.
The Exciton example has the more moderate dependence D = 3(2L + 1) 3 (see Table 3).Here, the parameter L is a geometric cutoff that limits the maximal distance between the electron and hole that constitute the exciton.This example has a number of other parameters that are adapted literally from [2].These parameters enter into the matrix entries, and thus affect the matrix spectrum and, finally, the algorithmic hardness of computing the eigenvalues of interest that determine the physical properties of the exciton.
Both Hubbard and Exciton are examples of difficult matrices, albeit for different reasons.For Hubbard, one unresolved challenge is to compute multiple interior eigenvalues for large n_fermions, n_sites, which becomes extremely difficult because of the rapid growth of the matrix dimension (specialized techniques for the Hubbard model such as the density-matrix renormalization group [70] cannot compute interior eigenvalues).Due to the irregular sparsity pattern of the Hubbard matrices (see Fig. 19 below), already the communication overhead of spMVM poses a serious obstacle to scalability and parallel efficiency.For Exciton, which are essentially stencil-like matrices of moderate size, the challenge is to compute some hundred eigenvalues out of a strongly clustered spectrum.Here, it is the  3 Matrix dimension D for the Hubbard and Exciton example, as a function of the respective parameter n_sites (and default value n_fermions = 5) or L Fig. 19 Sparsity pattern of the Hubbard (Hubbard,n_sites=40,n_fermions=20) and spin chain (SpinChainXXZ,n_sites=32,n_up=8) example poor convergence of iterative eigenvalue solvers for nearly degenerate eigenvalues that renders this problem hard.Thanks to the algorithmic advances in the ESSEX project, we now have reached a position that allows for future progress on these problems.
ScaMaC comes with several convenient features.For example, the Hubbard matrix includes the parameter ranpot to switch on a random potential.Random numbers in ScaMaC are entirely reproducible, and independent of the number of threads or processes that call the ScaMaC routines, or of the order in which the matrix rows are generated.An identical random seed gives the same matrix under all circumstances.In particular, individual matrix rows can be reconstructed at any time, which simplifies a fault-tolerant program design (see Sect. 5.3).Another feature is the possibility to effortlessly generate the (conjugate) transpose of nonsymmetric (non-Hermitian) matrices, which is considerably easier than constructing the transpose of a (distributed) sparse matrix after generation.

Eigensolvers in Quantum Physics: Graphene, Topological Insulators, and Beyond
Because of the linearity of the Schrödinger equation, quantum physics is a paradigm for numerical linear algebra applications.Historically, some application cases, such as the computation of the ground state (i.e. of the eigenvector to the minimal eigenvalue), have received so much attention that only gradual progress remains possible nowadays.In the ESSEX project we instead address two major cases where novel algorithmic improvements and systematic utilization of large-scale computing resources through state-of-the-art implementations still result in substantial qualitative progress.These two cases are the computation of (i) extreme eigenvalues with high degeneracy, which is addressed with a block Jacobi-Davidson algorithm, (ii) multiple interior eigenvalues, which is addressed by various filter diagonalization techniques.Application case (i) has been documented in [64], including the example of spin chain matrices (SpinChainXXZ in the ScaMaC).
For application case (ii) the primary quantum physics example are graphene [13] and topological insulators [31] (Graphene and TopIns in the ScaMaC).For these examples, eigenvalues towards the center of the spectrum, near the Fermi energy of the material, are those of interest.This situation is similar to applications in quantum chemistry and density functional theory, but in our case the matrices represent a full (disordered or structured) two or three-dimensional domain, and are usually larger than those considered elsewhere [57].
Starting with the paper [62] on Chebyshev filter diagonalization (ChebFD) and culminating in the BEAST software package (see Sect. 3.2), the computation of interior eigenvalues of large-to-huge graphene and topological insulator matrices has been successfully demonstrated with ESSEX algorithms, using polynomial filters derived from Chebyshev polynomials.Already with the simple ChebFD algorithm we could compute N T 100 eigenvectors from the center of the spectrum of a matrix with dimension D 10 9 (i.e. an effective problem size N T × D 10 11 ), in order to understand the electronic properties of a structured topological insulator (see Figure 13 in [62]).With improved filter coefficients and a more sophisticated implementation, the polynomial filters in the (P-) BEAST package deal with such problems at reduced computational cost (see Sect. 3.2).Such large-scale computations heavily rely on the optimized spMMVM and TSMM kernels of the GHOST library (see Sect. 5).
To appreciate the numerical progress reflected in these numbers one should note the different scaling of the numerical effort N MVM (measured in terms of the dominant operation of spMVM) for the computation of extreme and interior eigenvalues (cf. the discussion in [62]).In an idealized situation with equidistant eigenvalues, we have roughly N MVM ∼ D 1/2 for extreme but N MVM ∼ D for interior eigenvalues.For the D 10 9 example, we have to compensate for a factor 10 4 -10 5 to enable computation of interior instead of extreme eigenvalues.
Algorithm and software development in ESSEX has been to a large degree application-driven.Now, at the end of the ESSEX project, where the algorithms for our main application cases have become available, we follow two ways to go beyond the initial quantum physics applications.First, entirely new applications can now be addressed with ESSEX software, extending our efforts to non-linear and non-Hermitian problems (see Sect. 6.2).Second, relevant applications such as the Hubbard and Exciton examples (see Sect. 5.5) still fit into the two major application areas already addressed in ESSEX, but further increase the computational complexity.For Exciton, the strongly clustered spectrum with many nearly-degenerate eigenvalues leads to a numerical effort N MVM D 1/2 already for extreme eigenvalues.For Hubbard, the huge matrix dimension D is a serious obstacle for the computation of interior eigenvalues.
The Hubbard matrices also hint at an application-specific issue of general interest that we encountered but could not solve within ESSEX.Specifically, it is the complicated sparsity pattern of many of our quantum physics matrices (see Fig. 19) that adversely affects the parallel efficiency of distributed spMVM, and thus of our entire software solutions.Node-level performance engineering is here easily overcompensated by communication overhead.Unfortunately, the communication overhead is not reduced by standard matrix reordering strategies [61,71,84].This problem can be partially alleviated by overlapping communication with computation, as in the spMMVM (see Sect. 2), but a full solution to restore parallel efficiency is not yet available.Clearly, our different application scenarios still provide enough incentive to think about future numerical, algorithmic, and computational developments beyond the ESSEX project.

New Applications in Nonlinear Dynamical Systems
The block Jacobi-Davidson QR eigensolver in PHIST is capable of solving nonsymmetric and generalized eigenvalue problems of the form where B should be symmetric and positive definite.In [75], we exploited several unique features of this implementation to study the linear stability of a three-dimensional reaction-diffusion equation: the Jacobian is non-symmetric, the preconditioner was implemented in Epetra (which can be used directly as a backend for PHIST), and the high degree of symmetry in the model yields eigenvalues with high geometric multiplicity (up to 24).We therefore use a relatively large block size of 8 for these computations to achieve convergence to the desired 20-50 eigenpairs (λ i , x i ), with the real part of λ i near 0. In a recent Ph.D. thesis [74], the solver was also used for studying the linear stability of incompressible flow problems.Here B is in fact only semi-definite, and the preconditioner has to make sure that the solution stays in the 'divergence-free space', in which the velocity field satisfies ∇ • u = 0 and B induces a norm.Another ongoing effort concerning dynamical systems is the use of PHIST to parallelize the dynamical systems analysis tool PyNCT, which has as its main application the study of superconductors [82].We have taken first steps to use PHIST as backend for the Python-based algorithms in PyNCT.Furthermore, it is possible to solve the eigenvalue problems arising in PyNCT directly by the BJDQR method in PHIST.Our goal here is the scalable parallel and fully automatic computation of bifurcation diagrams using PyNCT and any backend supported by PHIST.
The Statistical Learning Lab led by Dr. Marina Meila at the University of Washington started to use the PHIST eigensolver to compute spectral gaps for Laplacian matrices obtained from conformation trajectories in molecular dynamics simulations, and other scientific data [14,59].These are symmetric positive definite matrices whose dimensions equal the number of simulation steps, typically of the order of n = 10 6 .When the data intrinsic dimension d is fixed, and much smaller than n (in our examples d < 10 ), the Laplacian is a sparse matrix.The sparsity pattern is not regular, and it is data dependent, as it reflects the neighborhood relationships in the data.Hence, in densely sampled regions rows will have many more non-zeros than in the sparsely populated regions of the data.In a manifold embedding algorithm, the eigengaps identify the optimal number of coordinates in which to embed the data.Furthermore, for data sizes n 10 6 , PHIST is used to compute the diffusion map embedding itself for the higher frequency coordinates for which existing methods are prohibitively slow.

International Collaborations
The internationalisation effort in the second phase of SPPEXA has fostered the ESSEX-II activities in several directions.First and foremost it amplified the scientific expertise in the project.Soon it became clear that complementing knowledge and developments could be leveraged across the partners.A specific benefit of the collaboration between German and Japanese partners is their very different background in terms of HPC infrastructures.Through close personal collaboration within the project all partners could easily access and use latest supercomputers on either side (see Sect. 4.2 and note that the BEAST framework has also been ported to the K-computer).Together with the joint collaboration on scientific problems and software development a steady exchange evolved with many personal research visits which also opened up collaboration with partners not involved directly in ESSEX-II.
The results described in Sect.3.1 on preconditioners are a direct result of the collaboration of ESSEX-II with the ppOpen-HPC 18 project led by Univ. of Tokyo.On the other hand, the CRAFT library developed at Univ. of Erlangen is utilized in an FEM code of Univ. of Tokyo and is part of a follow on JHPCN project with the German partner involved as associated partner.
Collaboration between Japanese and German working groups made possible the expansion of the BEAST framework for projection based eigensolvers to include Sakurai-Sugiura methods.Various numerical and theoretical issues associated with the implementation of the solver within an iterative framework were resolved, and new ideas explored during research visits.Results based on this collaboration have so far been presented in multiple conferences and a paper in preparation [35].
The linear systems arising from numerical quadrature in the BEAST-C and BEAST-M framework were used in the testing and development of a Block Cholesky-based ILU preconditioner.The integration of an interface to this solver into BEAST has begun.Examining strategies and expectations for solving these extreme ill-conditioned problems was a point of intense discussion and collaboration between working groups.One results was the development of RACE (see Sect. 3.4).Beyond the discussion between several Japanese and German ESSEX-II partners also a strong collaboration with the Swiss partner (O.Schenk) of EXASTEEL-II evolved, who is an expert on direct solvers and graph partitioning.In this context also a collaboration with T. Iwashita (Hokkaido Univ., Japan) started in terms of hardware efficient coloring.
Throughout the project, the variety of large matrices continuously added to the ScaMaC library allowed for testing with a variety of realistically challenging problems of both real and complex types in all ESSEX-II working groups.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/),which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Fig. 3 Fig. 4
Fig. 3 Computational time and convergence of BIC-COCG for a graphene benchmark problem (strong scaling)

Fig. 6
Fig. 6 Strong scaling of BEAST and z-Pares for a 1M × 1M standard eigenproblem based on a graphene sheet of dimension 2000 × 500.Both solvers found 260 eigenpairs in the interval [−0.01, 0.01] to a tolerance of 1 × 10 −8 .Both methods used 4 moments and began with random initial block vector Y .For BEAST, Y contained 100 columns; for z-Pares, 130. Testing performed on the Emmy HPC cluster at RRZE.MUMPS was used for the direct solution of all linear systems.N refers to the number of quadrature nodes along a circular contour, NP to the number of processes

Fig. 7
Fig.7Left: average residual over the BEAST iterations for using double or single precision throughout, and for switching from single to double precision in the 7th, 9th, or 11th iteration, respectively.Right: time (in seconds) to convergence for a size 1,048,576 topological insulator (complex values) with a search space size of 256 and a polynomial degree of 135 on 8 nodes of the Emmy-cluster at RRZE.Convergence is reached after identical numbers of iterations (with the exception of pure single precision, of course).The timings can vary for different ratios of polynomial degree and search space size and depend on the single precision performance of the underlying libraries

Fig. 8 A
Fig. 8 A map of the relationships among the contour integral-based eigensolvers

2 F 6 SFig. 9 Fig. 10
Fig.9 SymmSpMV performance of RACE compared to other methods.The roofline model for SymmSpMV is shown in Fig.9afor reference.Representative matrices from[17] and ScaMaC (see Sect. 5.5) were used.Note that the matrices are ordered according to increasing number of rows.(One Skylake Platinum 8160 CPU [24 threads]).(a) Performance of RACE compared with MKL.(b) Performance of RACE compared to other coloring approaches

Fig. 14
Fig.14 Weak scaling of BEAST-P on OFP, Piz Daint, and SNG, and strong scaling on SNG.Dashed lines denote ideal scaling with respect to the smallest number of nodes in the set.(a) Weak scaling of BEAST-P on OFP for problems of size 220 (2 nodes) to 230 (2048 nodes, about one quarter of the full machine).(b) Weak scaling of BEAST-P on Piz Daint for problems of size 221  (1 node) to 232 (2048 nodes).(From[56]).(c) Weak scaling of BEAST-P on SNG for problems of size 221 (1 node) to 1.53 × 232 (3136 nodes, about half of the full machine).(d) Strong scaling of BEAST-P on SNG for problems of size 2 28 (crosses) and 230 (triangles)

4 Fig. 15
Fig.15 Weak scaling behavior of the PHIST BJDQR solver for a symmetric PDE benchmark problem and different block sizes

Fig. 17
Fig.17 Lanczos application with various checkpoint/restart and process failure recovery scenarios using 128 nodes (256 processes) on the RRZE Emmy cluster.On average the communication recovery time is 2.6 s (ULFM-MPI v1.1)

Fig. 18
Fig. 18 Code example for row-by-row matrix generation with the generic ScaMaC generators

Table 1
The parameter values for Lanczos, JD, and Beast benchmarks

Table 2
Parameters of the Hubbard and Exciton matrix generator in the ScaMaC