Introduction to StarNEig -- A Task-based Library for Solving Nonsymmetric Eigenvalue Problems

TL;DR

StarNEig is a task-based library built on StarPU for solving dense nonsymmetric eigenvalue problems efficiently on shared and distributed memory systems, with GPU support planned.

Contribution

It introduces a new task-based library for nonsymmetric eigenvalue problems, offering compatibility with ScaLAPACK and supporting distributed and GPU-accelerated computations.

Findings

Demonstrates competitive performance in computational experiments.

Provides a ScaLAPACK compatibility layer for easy transition.

Supports both shared and distributed memory architectures.

Abstract

In this paper, we present the StarNEig library for solving dense non-symmetric (generalized) eigenvalue problems. The library is built on top of the StarPU runtime system and targets both shared and distributed memory machines. Some components of the library support GPUs. The library is currently in an early beta state and only real arithmetic is supported. Support for complex data types is planned for a future release. This paper is aimed for potential users of the library. We describe the design choices and capabilities of the library, and contrast them to existing software such as ScaLAPACK. StarNEig implements a ScaLAPACK compatibility layer that should make it easy for a new user to transition to StarNEig. We demonstrate the performance of the library with a small set of computational experiments.

Figures

Tables

Table 1. Table 1 : Current status of the StarNEig library.

Step	Shared memory	Distr. memory	GPUs
Hessenberg	Complete	ScaLAPACK	Single GPU
Schur	Complete	Complete	Experimental
Reordering	Complete	Complete	Experimental
Eigenvectors	Complete	In progress	Not planned
Hessenberg-triangular	Planned (LAPACK)	ScaLAPACK	Not planned
Generalized Schur	Complete	Complete	Experimental
Generalized reordering	Complete	Complete	Experimental
Generalized eigenvectors	Complete	In progress	Not planned

Table 2. Table 2 : A run time comparison between ScaLAPACK and StarNEig.

	CPU cores		Schur reduction (secs)		Eigenvalue reordering (secs)
$n$	ScaLAPACK	StarNEig	PDHSEQR	StarNEig	PDTRSEN	StarNEig
10 000	36	28	38	18	12	3
20 000	36	28	158	85	72	25
40 000	36	28	708	431	512	180
60 000	121	112	992	563	669	168
80 000	121	112	1667	904	1709	391
100 000	121	112	3319	1168	3285	737
120 000	256	252	3268	1111	2902	581

Equations

$$A{x}_i={\lambda }_i{x}_i,x_i\ne 0.$$

$$A{x}_i={\lambda }_{i}B{x}_i,x_i\ne 0.$$

$$A=Q_{1}H{Q}_1^T\text{ or }(A,B)=Q_1(H,R)Z_1^T,$$

$$H=Q_{2}S{Q}_2^T\text{ or }(H,R)=Q_2(S,T)Z_2^T,$$

$$(S-{\lambda }_{i}I)y_i=0\text{ or }(S-{\lambda }_{i}T)y_i=0$$

$$x_i=Q_1{Q}_2{y}_i\text{ or }{x}_i=Z_1{Z}_2{y}_i.$$

$$S=Q_3\hat{S}{Q}_3^T\text{ or }(S,T)=Q_3(\hat{S},\hat{T})Z_3^T,$$

Full text

1 Introduction

In this paper, we present the StarNEig library [1] for solving dense non-symmetric (generalized) eigenvalue problems. StarNEig differs from the existing libraries such as LAPACK and ScaLAPACK in that it is relies on a modern task-based approach. More specifically, StarNEig is build on top of the StarPU runtime system [2]. This allows StarNEig to target both shared memory and distributed memory machines. Furthermore, some components of the StarNEig library support GPUs. The library is currently in an early beta state and under under continuous development.

This paper targets potential users of the library. We hope that readers, who are already familiar with ScaLAPACK, will be able to decide if StarNEig is suitable for them. In particular, we want to communicate what type of changes are necessary to make their software work with StarNEig. We will explain, through an example, why the task-based approach can potentially lead to superior performance when compared to older, but well-established, approaches. We also present a small sample of computational results which demonstrate the expected performance of the library. We refer the reader to [13] for more comprehensive performance and accuracy evaluations.

The rest of this paper is organized as follows: Section 2 provides a brief introduction to the solution of dense eigenvalue problems. Section 3 explains the task-based approach and Section 4 introduces the reader to some of the inner workings of StarNEig. Section 5 presents a small set of computational results and, finally, Section 6 concludes the paper.

2 Solution of Dense Nonsymmetric Eigenvalue Problems

Given a matrix $A\in\mathbb{R}^{n\times n}$, the standard eigenvalue problem consists of computing eigenvalues $\lambda_{i}\in\mathbb{C}$ and matching eigenvectors $x_{i}\in\mathbb{C}^{n}$ such that

[TABLE]

Similarly, given matrices $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times n}$ the generalized eigenvalue problem for the matrix pair $(A,B)$ consists of computing generalized eigenvalues $\lambda_{i}\in\mathbb{C}$ and matching generalized eigenvectors $x_{i}\in\mathbb{C}^{n}$ such that

[TABLE]

If the matrix $A$ or the matrices $A$ and $B$ are dense and nonsymmetric, then route of acquiring the (generalized) eigenvalues and the (generalized) eigenvectors usually includes the following three steps:

Hessenberg(-triangular) reduction:

The matrix $A$ or the matrix pair $(A,B)$ is reduced to upper Hessenberg or Hessenberg-triangular form by an orthogonal similarity transformation

[TABLE]

where $H$ is upper Hessenberg, $R$ is a upper triangular, and $Q_{1}$ and $Z_{1}$ are orthogonal.

Schur reduction:

The Hessenberg matrix $H$ or the Hessenberg-triangular matrix pair $(H,R)$ is reduced to Schur or generalized Schur form by an orthogonal similarity transformation

[TABLE]

where $S$ is upper quasi-triangular with $1\times 1$ and $2\times 2$ blocks on the diagonal, $T$ is a upper triangular, and $Q_{2}$ and $Z_{2}$ are orthogonal. The eigenvalues or generalized eigenvalues can be determined from the diagonal blocks of $S$ or $(S,T)$.

Eigenvectors:

Finally, we compute vectors $y_{i}\in\mathbb{C}^{n}$ from

[TABLE]

and backtransform to the original basis by

[TABLE]

Additionally, a fourth step can be performed to acquire a desired invariant subspace of $A$ or $(A,B)$:

Eigenvalue reordering:

The Schur form $S$ or the generalized Schur form $(S,T)$ is reordered, such that a selected set of eigenvalues or generalized eigenvalues appears in the leading diagonal blocks of an updated Schur form $\hat{S}$ or an updated generalized Schur form $(\hat{S},\hat{T})$, by an orthogonal similarity transformation

[TABLE]

where $Q_{3}$ and $Z_{3}$ are orthogonal.

See [4] for a detailed explanation of the underlying mathematical theory.

3 A Case for the Task-Based Approach

A task-based algorithm functions by cutting the computational work into self-contained tasks that all have a well defined set of inputs and outputs. The algorithm inserts the tasks into a runtime system that derives the task dependences and schedules the tasks to computational resources in a sequentially consistent order. As long as the cutting is carefully done, the underlying parallelism is exposed automatically as the runtime system unravels the resulting task graph.

3.1 Novelty in StarNEig

The first main source of novelty in StarNEig comes from the way in which the computational work is cut into tasks. As with many other task-based matrix algorithms, the matrices are divided into (square) tiles and each task takes a set of tiles as its input and produces/modifies a set of tiles as its output. The Hessenberg reduction, Schur reduction and eigenvalue reordering steps are based on two-sided transformation algorithms. These algorithms lead to data dependences that are much more complicated than the dependences arising from one-sided transformation algorithm such as the LU factorization. The second main source of novelty in StarNEig is related to the eigenvector computation step. Here, the data dependences are comparatively simple but the computations must be protected against floating-point overflow. This is a nontrivial issue to address in a parallel setting; see [7, 8, 9].

Furthermore, the Schur reduction and eigenvalue reordering steps apply a series of overlapping local transformations to the matrices. Due to this overlap, the two computational steps cannot have have a clear one-to-one mapping between the tasks and the (output) tiles since the local transformations must at some point cross between two or more tiles. Instead, most task ends up modifying several tiles and this can introduce spurious data dependences111A spurious data dependency is created when two (or more) tasks modify non-overlapping parts of the same tile but the runtime system interprets this as a true data dependency. that limit the concurrency.

3.2 Bulge Chasing and Eigenvalue Reordering

We will now use the Schur reduction and eigenvalue reordering steps to illustrate some benefits of the task-based approach. The modern approach for obtaining a Schur form $S$ of $A$ is to apply the multishift QR algorithm with Aggressive Early Deflation (AED) to the upper Hessenberg form $H$ (see [6] and references therein). The algorithm is a sequence of steps of two types: AED and bulge chasing. The bulge chasing step creates a set of bulges which are chased down the diagonal to complete one pipelined QR iteration. This is accomplished by applying sequences of overlapping $3\times 3$ Householder reflectors to $H$. Similarly, the eigenvalue reordering step is based on applying sequences of overlapping Givens rotations and $3\times 3$ Householder reflectors to $S$.

If the local transformations are applied one by one, then memory is accessed as shown in Fig. 1(a). This is grossly inefficient for two reasons: i) the transformation is so localized that parallelizing it would not produce any speedup and ii) the matrix elements are touched only once thus leading to very low arithmetic intensity. The modern approach groups a set of local transformation together and initially applies them to a relatively small diagonal window as shown in Fig. 1(b). The local transformations are accumulated into an accumulator matrix and later applied as level-3 BLAS operations acting on the relevant sections of the matrix. This leads to much higher arithmetic intensity and enables proper parallel implementations as multiple diagonal windows can be processed simultanously.

In particular, the Schur reduction and eigenvalue reordering steps are implemented in ScaLAPACK as PDHSEQR [6] and PDTRSEN [5] subroutines, receptively. Following a typical ScaLAPACK formula, the matrices are distributed in two-dimensional block cyclic fashion. The resulting memory access pattern is illustrated in Fig. 1(c) for a $3\times 3$ MPI process mesh. In this particular example, three diagonal windows can be processed simultaneously. The related level-3 BLAS updates require careful coordination since the left and right hand side updates must be performed in a sequentially consistent order. In practice, this means (global or broadcast) synchronization after each set of (left or right hand side) updates have been applied.

In a task-based approach, this can be done using the following task types:

Window task

applies a set of local transformations inside the diagonal window. Takes the intersecting tiles as input, and produces updated tiles and an accumulator matrix as output.

Right update task

applies accumulated right-hand side updates using level-3 BLAS operations. Takes the intersecting tiles and an accumulator matrix as input, and produces updated tiles as output.

Left update task

applies accumulated left-hand side updates using level-3 BLAS operations. Takes the intersecting tiles and an accumulator matrix as input, and produces updated tiles as output.

The tasks are inserted in a sequentially consistent order and each window chain leads to a task graph like the one shown in Fig. 2. It is critical to realize that the runtime system guarantees that the tasks are executed in a sequentially consistent order. In particular, there is no need for synchronization and different stages are allowed to overlap as illustrated in Fig. 1(d). This leads to much higher concurrency. Under suitable conditions, the AED step can also be partially overlapped with the bulge chasing step. Other benefits of the task-based approach include, for example, better load balancing, task priorities, accelerators support and implicit MPI communication. See [10, 11, 12, 13] for implementation details and further information.

4 StarNEig Library

StarNEig is a C-library that runs on top of the StarPU task-based runtime system. StarPU handles low-level operations such as heterogeneous scheduling, data transfers and replication between various memory spaces, and MPI communication between compute nodes. In particular, StarPU is responsible for managing the various computational resources such as CPU cores and GPUs. In order to accomplish this, StarPU creates a set of worker threads; usually one thread per computational resource. In addition, one thread is responsible for inserting the tasks and tracking the state of the machine. If necessary, one additional thread is allocated for MPI communication. Thus, StarNEig should be used in a one process per node (1ppn) configuration, i.e., several CPU cores should be allocated for each process (a node can be a full compute node, a NUMA island or some other collection of CPU cores).

The current status of the StarNEig library is summarized in Table 1. The library is currently in an early beta state. At the time of writing this paper, only real arithmetic is supported and certain interface functions are implemented as LAPACK and ScaLAPACK wrapper functions. The library is under continuous development. In particular, additional distributed memory functionality and support for complex data types are planned for a future release.

4.1 Distributed Memory

StarNEig distributes the matrices in rectangular blocks of a uniform size (excluding the last block row and column) as illustrated in Fig. 3(a). The data distribution, i.e., the mapping from the distributed blocks to the MPI process rank space, can be arbitrary as illustrated in Fig. 3(b). A user has three options:

Use the default data distribution. This is recommended for most users and leads to reasonable performance in many cases. 2. 2.

Use a two-dimensional block cyclic distribution (see Fig. 3(c)). In this case, the user may select the MPI process mesh dimensions and the rank ordering. 3. 3.

Define a data distribution function $d:\mathbb{Z^{+}}\times\mathbb{Z^{+}}\to\mathbb{Z^{+}}$ that maps the block row and column indices to the MPI rank space. For example, in Fig. 3(b), the rank 0 owns the blocks (0,1), (1,2), (1,5), (1,6), (2,6), (3,0) and (3,5).

The library implements distribution agnostic copy, scatter and gather operations.

Users who are familiar with ScaLAPACK are likely accustomed to using relatively small distributed block sizes (between 64–256). In contrast, StarNEig functions optimally only if the distributed blocks are relatively large (at least 1000). This is due to the fact that StarNEig further divides the distributed blocks into tiles and a tiny tile size leads to excessive task scheduling overhead because the tile size is closely connected to the task granularity. Furthermore, as mentioned in the preceding section, StarNEig should be used in 1ppn configuration as opposed to a one process per core (1ppc) configuration which is more common with ScaLAPACK.

4.2 ScaLAPACK Compatibility

StarNEig is fully compatible with ScaLAPACK and provides a ScaLAPACK compatibility layer that encapsulates BLACS contexts and descriptors [3] inside transparent objects, and implements a set of bidirectional conversion functions. The conversions are performed in-place and do not modify any of the underlying data structures. Thus, users can mix StarNEig interface functions with ScaLAPACK subroutines without intermediate conversions.

5 Performance Evaluation