AICS TECHNICAL REPORT
N0. 2015-002

USER’S MANUALEIGEN EXAVERSION 2.3C
By

EIGEN EXA DEVELOPMENT GROUPLARGE-SCALE PARALLEL NUMERICAL COMPUTING TECHNOLOGY RESEARCH TEAMRIKEN ADVANCED INSTITUTE FOR COMPUTATIONAL SCIENCE

Submitted on 28/08/2015
Accepted on 04/09/2015
Published and copyrighted by
RIKEN Advanced Institute for Computational Science (AICS)
7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, 650-0047, Japan

User’s Manual

EigenExa
Version 2.3c
EigenExa Development Group
Large-scale Parallel Numerical Computing Technology Research TeamRIKEN Advanced Institute for Computational Science
24 June 2015

2

Contents
1 Introduction
1.1 EigenExa and its course of development . . . . . . . . . . . . . . . . . . . . . . .
1.2 License for use and copyright . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55
62 Before use
2.1 Software required for EigenExa installation . . . . . . . . . . . . . . . . . . . . .
2.2 Obtaining EigenExa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Compile and install procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
7
7

3 Quick tutorial

9

4 API
13
4.1 eigen init . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.2 eigen free . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.3 eigen get blacs context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.4 eigen sx . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.5 eigen s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.6 eigen get version . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.7 eigen show version . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.8 eigen get matdims . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.9 eigen memory internal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.10 eigen get comm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.11 eigen get procs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.12 eigen get id . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.13 eigen loop start . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.14 eigen loop end . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.15 eigen translate l2g . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.16 eigen translate g2l . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.17 eigen owner node . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.18 KMATH EIGEN GEV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
5 Other key considerations
19
5.1 Caution regarding compatibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.2 Binding with other languages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.3 Behavior on error occurrence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.4 Shared library handling in versions 1.x . . . . . . . . . . . . . . . . . . . . . . . . 19
3

4

CONTENTS

A Algorithm overview
21
A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
A.2 Various approaches and related projects . . . . . . . . . . . . . . . . . . . . . . . 21
A.3 eigen s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
A.4 eigen sx . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
A.5 Diﬀerences between eigen s and eigen sx . . . . . . . . . . . . . . . . . . . . . . . 24
A.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
Acknowledgements

27

References

29Chapter 1

Introduction
1.1

EigenExa and its course of development

EigenExa is a high-performance numerical eigenvalue solver. It traces its history back to EigenES(code name; no formal name) [1] developed on the Earth Simulator, which attained the rank of
number one in the world of such systems in 2002. EigenES was nominated for a Gordon Bell Prizeat SC 2006 and today continues to serve as an eigenvalue solver on large-scale PC clusters. This
led to the initiation of EigenK [2, 3] development around 2008. The EigenK library became theimmediate predecessor of EigenExa, and in August 2013 EigenK was renamed EigenExa andpublic release was begun, following entry into operation on the K computer[4, 5]. EigenExa
development continues, with the underlying objective being to achieve an eigenvalue libraryscalable to operate on future post-petascale (“exa” or “extreme”) computer systems.
In its present release (version 2.3c), EigenExa provides the simplest function of computing
all eigenpairs (eigenvalues paired with their respective eigenvectors) for both standard and generalized eigenvalue problems. As reported elsewhere [2, 3], EigenExa applies both classical andadvanced algorithms in the same basic manner as EigenK, and thereby reduces the required
computation time for diagonalization.
The development of EigenExa includes the utilization of various parallel programming languages and libraries, encompassing MPI, OpenMP, high-performance BLAS, and SIMD vectorized Fortran90 compiler techniques. EigenExa is expected to open the way for high-performancecomputation through multiple simultaneous functions characterized by the following.
1. Inter-node parallelism in distributed memory architecture, by MPI
2. Parallelism in shared-memory parallel computers and multicore processors, by OpenMP
3. High parallelism utilizing BLAS highly optimized by vendors
4. SIMD or coarse-grained parallelism utilizing vendor-provided high-performance compilers
The good features of Fortran90 are also actively incorporated into EigenExa. The API of
EigenExa is more ﬂexible than libraries implemented in Fortran77, and it provides a user-friendly
interface, based on modular interfaces and optional parameters. Data distribution is limited to
two-dimensional cyclic division, the processor map can be speciﬁed in almost any arbitrary
conﬁguration, and compatibility and consistency with existing numerical computation libraries
are guaranteed if the data redistribution function provided by ScaLAPACK is used. EigenExa
also permits user interface speciﬁcation (or omission) for heightened performance, such as blockparameters that strongly aﬀect execution performance.
5

6

CHAPTER 1. INTRODUCTION

In the parallel performance of the library itself, it achieves heightened performance by reducing the EigenK communication overhead, and it has been shown that in most cases EigenExa
performance exceeds that of EigenK, ScaLAPACK, and others of the highest-level numericalcomputation libraries [3].
Today, EigenExa is in operation on many HPC platforms, including the K computer and
its Fujitsu PRIMEHPC FX10 commercial variant, various cluster computers using Intel x86
series processors, IBM Blue/Gene Q systems, and the NEC vector computer SX series systems.
Furthermore, reports on EigenExa have been presented at scientiﬁc conferences [6, 7, 8, 9], so ifnecessary, please refer to them.
This user’s manual for EigenExa version 2.3c covers procedures from installation to actual
use, with particular consideration given to installation and compiling, a quick tutorial, the API
list, and EigenK compatibility. It is written and provided with the hope of all EigenExa teamdevelopers that it will be helpful for achieving eﬃcient parallel simulations by many users.

1.2

License for use and copyright

Permission to use KMATH RANDOM is granted on the basis of the BSD 2-Clause License(found in LICENCE.txt in the library).

LICENCE.txt



Copyright (C) 2012- 2014 RIKEN.
Copyright (C) 2011- 2012 Toshiyuki Imamura
Graduate School of Informatics and Engineering,
The University of Electro-Communications.
Copyright (C) 2011- 2014 Japan Atomic Energy Agency.
-----------------------------------------------------------------------Copyright notice is from here-----------------------------------------------------------------------Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:* Redistributions of source code must retain the above copyrightnotice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyrightnotice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS‘‘AS IS’’ AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHTHOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
------------------------------------------------------------------------

Chapter 2

Before use
2.1

Software required for EigenExa installation

Several software packages are needed to compile the EigenExa library. BLAS, LAPACK, ScaLAPACK, and MPI must be installed in the system before EigenExa is compiled. To present, it
has been conﬁrmed that EigenExa can be compiled with the following libraries.
BLAS
LAPACK
ScaLAPACK
MPI

2.2Intel MKL, GotoBLAS, OpenBLAS, ATLAS
Fujitsu SSL II, IBM ESSL, NEC MathKeisan
Version 3.4.0 or later
Version 1.8.0 or later
MPICH2 version 1.5 or later, MPICH version 3.0.2 or later
OpenMPI version 1.6.4 or later, MPI/SX

Obtaining EigenExa

All available information on EigenExa can be obtained at the following URL.
http://www.aics.riken.jp/labs/lpnctrt/EigenExa.html
Tarball distribution is also performed via this URL. Planning is in progress for provision offurther information on EigenExa.

2.3

Compile and install procedure

Several steps are necessary to compile the EigenExa library. Proceed as described in the following
installation guideline.
Decompression and extraction First, unpack the tarball on the working directory and thenmove to the EigenExa-2.3c directory% tar zxvf EigenExa-2.3c.tgz% cd EigenExa-2.3c
7

8

CHAPTER 2. BEFORE USE

Environment setting Second, edit Makefile and make_inc.xxx to match the user’s environment. For xxx, select and enter the character string corresponding to the compiler to be
used, from among the following:
1. BX900
2. Intel or Intel.shared
3. K_FX10 or K_FX10.shared
4. gcc5. BlueGeneQ
6. SX
The “.shared” suﬃx is speciﬁed in creating a shared library.
make Third, execute make. This results in creation of the static library libEigenExa.a or
the shared library libEigenExa.so.
% makeinstall Lastly, copy the library itself, libEigenExa.a (or libEigenExa.so for the sharedlibrary), eigen_libs.mod, and eigen_blacs.mod to the install directory, and end the procedure.
To install /usr/local/lib, for example,% cp libEigenExa.a eigen_libs.mod eigen_blacs.mod/usr/local/lib/

generalized eigenvalue computation driver routine When the generalized eigenvalue
driver routine KMATH EIGEN GEV with the same number as the installed version is downloaded fromhttp://www.aics.riken.jp/labs/lpnctrt/KMATH_EIGEN_GEV.html
and followed by make, the KMATH_EIGEN_GEV.o driver module for generalized eigenvalues is created. Linking this object when program linking enables computation of generalized eigenvalues.
At present, this function is a higher-level independent module used as the EigenExa eigenvalue
computation engine. In this manual, it is described in terms of its relation to EigenExa.

Chapter 3

Quick tutorial
The standard benchmark code can be obtained by moving to the working directory and executing‘make benchmark’. The source code components ‘main2.F’ and ‘Makefile’ should be useful forcode creation. The kernel of main2.f is as follows.

main2.f

use MPIuse eigen_libs...
call MPI_Init_thread( MPI_THREAD_MULTIPLE, i, ierr )call eigen_init( )N=10000; mtype=0
call eigen_get_matdims( N, nm, ny )allocate ( A(nm,ny), Z(nm,ny), w(N) )call mat_set( N, a, nm, mtype )call eigen_sx( N, N, a, nm, w, z, nm, m_forward=32, m_backward=128)
deallocate ( A, Z, w )...
call eigen_free( )call MPI_Finalize( ierr )
end





The above code represents only the framework and does not actually operate, but is suﬃcient
to indicate the overall ﬂow from initialization to sequence assurance to eigenvalue computation
to termination procedure.
In the above example, the initialization function eigen_init() is executed by a parameteromitted call. In eigen_init(), the group that executes the eigenvalue computation can be
speciﬁed as a communicator by the form comm=XXX. For parallel execution of eigenvalue computation simultaneously by multiple groups, parallel computation is enabled by passing thecommunicator created by MPI_Comm_split(), for example. Note that a constraint exists in thecurrent implementation such that performance of collective operation with respect to the communicator MPI_COMM_WORLD within eigen_init() requires that eigen_init() be called by allprocesses belonging to MPI_COMM_WORLD simultaneously. Since a diﬀerent communicator can be
speciﬁed for each individual process, a process not participating in the eigenvalue computation
can therefore have MPI_COMM_NULL speciﬁed for eigen_init(), thus having the call skipped
for the eigenvalue computation driver eigen_sx() itself. In short, it is possible to perform
910

CHAPTER 3. QUICK TUTORIAL

simultaneous execution of various operations other than eigen_sx() processing that include
eigen_sx().


MPI Comm split and MPI COMM NULLif ( my_rank < 10 ) thencomm = MPI_COMM_SELF
elsecomm = MPI_COMM_NULL
endif
call eigen_init( comm )call eigen_sx( .... )In EigenExa, processes belonging to a communicator speciﬁed by eigen_init() are deployed
and used in a two-dimensional process grid. EigenExa is designed to utilize a process grid that isnearly square in shape, to reduce communication traﬃc to the minimum possible. To heighten
user convenience, EigenExa has been developed so that the two-dimensional Cartesian grid
adopted by MPI can be speciﬁed to comm. If the Cartesian shape is two dimensional, thenin principle EigenExa can be called and computation performed for any process array, thus
enabling execution of complicated parallel processing through combination of the above multitype communicators. Because the Cartesian process grid is essentially row-major, the Cartesianis prioritized in the event of conﬂict with an order=’C’ speciﬁcation. For historical reasons, thedefault process grid in EigenExa is column-major.
The generation of matrix data is performed in mat_set(), which is called just before
eigen_sx(). The matrix data are distributed on the speciﬁed two-dimensional process gridin two-dimensional cyclic division style, and are stored in each process as a local array. Because
only some of the data are stored for each process, when matrix elements are accessed, a rule fortransformation between global and local indices is required.
The following program is an excerpt from mat_set() and is presented to compare the program for generation of a Frank matrix with a global counterloop structure and the same buttranslated to local counterloops.

matset(before)! Global loop program to compute a Frank matrix
do i = 1, n
do j = 1, na(j, i) = (n+1-Max(n+1-i,n+1-j))*1.0D+00
end doend do↓↓↓↓↓↓

11
matset(after)! Translated local loop program to compute a Frank matrix
use MPIuse eigen_libs
call eigen_get_procs( nnod, x_nnod, y_nnod )call eigen_get_id( inod, x_inod, y_inod )j_2 = eigen_loop_start( 1, x_nnod, x_inod )j_3 = eigen_loop_end ( n, x_nnod, x_inod )i_2 = eigen_loop_start( 1, y_nnod, y_inod )i_3 = eigen_loop_end ( n, y_nnod, y_inod )do i_1 = i_2, i_3
i = eigen_translate_l2g( i_1, y_nnod, y_inod )do j_1 = j_2, j_3
j = eigen_translate_l2g( j_1, x_nnod, x_inod )a(j_1, i_1) = (n+1-Max(n+1-i,n+1-j))*1.0D+00
end doend doThe eigen_loop_start() and eigen_loop_end() are used to transform the loop range.
The second and third parameters specify the process number and process ID derived from thecommunicator, which shows the direction of distribution. In this manual, the correspondence
is always [row] → “x” and [column] → “y” (in the case of the communicator for the overall
participating process, the “x” and “y” portions are characterless). It is important to note thatin EigenExa process IDs are managed as integers, starting with 1. The process ID obtained by
the query function eigen_get_id( ) therefore diﬀers from the MPI rank by 1, and the ID must
be reduce by 1 in cases where the MPI rank is required.
In the above program, the local loop counter value is translated to the corresponding global
counter value to be used, with eigen_translate_l2g() used for the this translation. The second and third parameters should be speciﬁed in a manner similar to eigen_loop_start(),for example. Conversely, to convert the global counter value to the local counter value,
eigen_translate_g2l() is used, with the proviso that if the global counter value is viewed
as a loop value, then eigen_translate_g2l() returns the corresponding local counter valueon the process that becomes the owner process (the process where the local counter value corresponding with the global counter value is included in the loop). In order to determine theowner process of the speciﬁed global loop counter, eigen_owner_node()is used. It should be
used when broadcasting in reference to a particular row or column vector value.
Furthermore, when progressing to computation together with ScaLAPACK for personnel atadvanced levels, the process grid context used by EigenExa should be obtained via the auxiliary
function eigen_get_blacs_context(), referring to the mtype=2 portion of the mat_set() function (the following shows the kernel of the PDTRAN() call that stores the matrix AS transposein matrix A).

12

CHAPTER 3. QUICK TUTORIAL

pdtran! Cooperation with ScaLAPACKNPROW = x_nnod; NPCOL = y_nnod
ICTXT = eigen_get_blacs_context( )CALL DESCINIT( DESCA, n, n, 1, 1, 0, 0, ICTXT, nm, INFO )! A <-- AS^T
CALL PDTRAN( n, n, 1D0, as, 1, 1, DESCA, 1D0, a, 1, 1, DESCA )When compiling, in addition to the use of mpif90, it is necessary to set the path (in mostcases, the -I option) because of the need to access eigen_libs and other modules. To link the
EigenExa library, it is also necessary to simultaneously link MPI, OpenMP, ScaLAPACK (ifversion 1.8 or earlier, then also BLACS), and so forth. In the case of Intel-compiler-based MPI,
the procedure is as follows (note that the library names around ScaLAPACK and BLAS varywith the environment).
% mpif90 -c a.f -openmp -I/usr/local/include -I/usr/local/lib% mpif90 -o exe a.o -openmp -L/usr/local/lib -lEigenExa -lscalapack \-llapack -lblas

Chapter 4

API
This section lists the functions in ‘eigen_libs.mod’ that have been assigned a public attribute.
The ﬁrst three routines are the main drivers and the others are utility functions. Parameters
having optional attributes attached (written in italics) can be omitted, and can also be speciﬁedby TERM=variable or constant value in the Fortran format form.

4.1

eigen init

Initializes the functions of EigenExa. Process grid mapping can be speciﬁed via the parameter‘comm’ or ‘order’. This procedure is collective (because collective operation is performed for
MPI_COMM_WORLD internally in the present version, after which all processes must call this procedure); comm can specify a diﬀerent value for each process group, and when diﬀerent processgroups simultaneously call driver functions (eigen_sx() or eigen_s()), parallel operations are
performed in driver function units. If comm is MPI_COMM_NULL, calling the handlers eigen_sx()
and eigen_s() results in an immediate return with no internal action. An inter-communicator
cannot be used for comm.
subroutine eigen init( comm, order )1. integer, optional, intent(IN) :: comm = MPI_COMM_WORLD
Base communicator.
Process grid mapping is enabled when a two-dimensional Cartesian grid is
speciﬁed for comm.
Note: If omitted, then MPI COMM WORLD
2. character*(*), optional, intent(IN) :: order = ’C’
Row or Column
Note: If omitted, then treated as ‘C’. If the grid major conﬂicts
with the Cartesian comm speciﬁcation, then ‘R’ is used.

4.2

eigen freeTerminates EigenExa functions.
subsroutine eigen free( flag )1. integer, optional, intent(IN) :: flag = 0
Timer printer ﬂag.
This parameter is for development, and therefore not ordinarily speciﬁed.
If omitted, then 0.
13

14

4.3

CHAPTER 4. API

eigen get blacs context

Returns the ScaLAPACK (BLACS) context corresponding to the process grid information speciﬁed by EigenExa.
integer function eigen get blacs context( )4.4eigen sxThis is the main EigenExa driver routine. Eigenpairs are computed via transformation to apentadiagonal matrix. This is a collective operation driver; all processes attributed to theprocess group calling must participate in the call.
subroutine eigen sx( n, nvec, a, lda, w, z, ldz, \m forward, m backward, mode )1.
integer, intent(IN) :: nMatrix and vector dimensions2.
integer, intent(IN) :: nvec
Number of computed eigenvectors
At present, this option is not supported;eigen_sx() computes all eigenvectors.
3.
real(8), intent(INOUT) :: a(lda,*)Symmetric matrix to be diagonalized
Array content is destroyed upon subroutine termination but
the FLOPS count is stored in a(1, 1).
4.
integer, intent(IN) :: lda
Leading dimension of array a5.
real(8), intent(OUT) :: w(n)Eigenvalues in ascending order6.
real(8), intent(OUT) :: z(ldz,*)
Orthogonal eigenvectors of matrix a
7.
integer, intent(IN) :: ldz
Leading dimension of array z
8.
integer, optional, intent(IN) :: m_forward = 48
Householder transformation blocksize (must be even number);
48 when omitted
9.
integer, optional, intent(IN) :: m_backward = 128
Householder back transformation blocksize; 128 when omitted
10. character, optional, intent(IN) :: mode = ’A’‘A’ : all eigenvalues and corresponding eigenvectors (default)‘N’ : eigenvalues only
‘X’ : add to mode A to improve eigenvalue precision

4.5

eigen sThis is the EigenExa driver routine.

4.6. EIGEN GET VERSION

15subroutine eigen s( n, nvec, a, lda, w, z, ldz, \m forward, m backward, mode )1.
integer, intent(IN) :: nMatrix and vector dimensions2.
integer, intent(IN) :: nvec
Number of computed eigenvectors
At present, this option is not supported;eigen_s() computes all eigenvectors.
3.
real(8), intent(INOUT) :: a(lda,*)Symmetric matrix to be diagonalized
Array content is destroyed upon subroutine termination but
the FLOPS count is stored in a(1, 1).
4.
integer, intent(IN) :: lda
Leading dimension of array a5.
real(8), intent(OUT) :: w(n)Eigenvalues in ascending order6.
real(8), intent(OUT) :: z(ldz,*)
Orthogonal eigenvectors of matrix a
7.
integer, intent(IN) :: ldz
Leading dimension of array z
8.
integer, optional, intent(IN) :: m_forward = 48
Householder transformation blocksize; 48 when omitted
9.
integer, optional, intent(IN) :: m_backward = 128
Householder back transformation blocksize; 128 when omitted
10. character, optional, intent(IN) :: mode = ’A’‘A’ : all eigenvalues and corresponding eigenvectors (default)‘N’ : eigenvalues only
‘X’ : add to mode A to improve eigenvalue precision

4.6

eigen get version

Returns the version information of EigenExa.
subroutine eigen get version( version, data, vcode )1. integer, intent(OUT) :: version
Version number in three digits.
The digits indicate major version, minor version, and patch level from the highest.
2. character, intent(OUT) :: dateRelease date.
3. character, intent(OUT) :: vcode
Code name of the version.

4.7

eigen show version

Print the version information of EigenExa to standard output.
subroutine eigen show version()

4.8

eigen get matdims

Returns the array size recommended by EigenExa. It is desirable for the user to dynamicallyallocate the local array using the array dimensions (nx,ny) obtained by this function or a larger

16

CHAPTER 4. API

numeric value. The overall matrix is (CYCLIC, CYCLIC) divided.
subroutine eigen get matdims( n, nx, ny )1. integer, intent(IN) :: nMatrix dimensions
2. integer, intent(OUT) :: nx
Lower limit of leading dimensions of arrays a and z
3. integer, intent(OUT) :: ny
Lower limit of second indices of arrays a and z

4.9

eigen memory internalThis function returns the dynamically allocated internal memory size, while EigenExa is called.
The user should ascertain the value returned by this function to ensure avoidance of memory
shortage. The type of the return value is changed to integer(8) in version 2.3c. If the return
value is −1 (a negative value), this funtion warns of possible integer overﬂow resulting from toolarge matrix dimension in EigenExa
integer(8) function eigen memory internal( n, lda, ldz, m1, m0 )1. integer, intent(IN) :: nMatrix dimensions
2. integer, intent(IN) :: lda
Leading dimension of array a
3. integer, intent(IN) :: ldz
Leading dimension of array z
4. integer, intent(IN) :: m1
Householder transformation blocksize (must be even number)
5. integer, intent(IN) :: m0
Householder transformation blocksize

4.10

eigen get comm

Returns the MPI communicator generated by eigen_init().
subroutine eigen get comm( comm, x comm, y comm )1. integer, intent(OUT) :: commBase communicator.
2. integer, intent(OUT) :: x_comm
Row communicator, attribute of all processes with matching row id.
3. integer, intent(OUT) :: y_comm
Column communicator, attribute of all processes with matching column id.

4.11

eigen get procs

Returns information on the number of processes related to the communicator generated byeigen_init().
subroutine eigen get procs( procs, x procs, y procs )1. integer, intent(OUT) :: procs
Number of processes in comm2. integer, intent(OUT) :: x_procs
Number of processes in x_comm
3. integer, intent(OUT) :: y_procs
Number of processes in y_comm4.12. EIGEN GET ID

4.12

17

eigen get id

Returns information on the process ID related to the communicator generated by eigen_init().
Here, the process ID diﬀers from the MPI rank; it is an integer value starting from 1 in therelation MPI rank = process ID − 1.
subroutine eigen get id( id, x id, y id )1. integer, intent(OUT) :: id
Process ID deﬁned by comm2. integer, intent(OUT) :: x_id
Process ID deﬁned by x_comm
3. integer, intent(OUT) :: y_id
Process ID deﬁned by y_comm

4.13

eigen loop start

Returns the loop starting value in the local loop structure corresponding to the speciﬁed global
loop starting value.
integer function eigen loop start( istart, nnod, inod )1. integer, intent(IN) :: istart
Global loop starting value
2. integer, intent(IN) :: nnod
Process number
3. integer, intent(IN) :: inod
Process ID

4.14

eigen loop end

Returns the terminal value of the loop in the local loop structure corresponding to the terminal
value of the speciﬁed global loop.
integer function eigen loop end( iend, nnod, inod )1. integer, intent(IN) :: istart
Global loop terminal value
2. integer, intent(IN) :: nnod
Process number
3. integer, intent(IN) :: inod
Process ID

4.15

eigen translate l2g
integer function eigen translate l2g( ictr, nnod, inod )1. integer, intent(IN) :: ictr
Local counter
2. integer, intent(IN) :: nnod
Process number
3. integer, intent(IN):: inod
Process ID

18

4.16

CHAPTER 4. API

eigen translate g2l
integer function eigen translate g2l( ictr, nnod, inod )1. integer, intent(IN) :: ictr
Global counter
2. integer, intent(IN) :: nnod
Process number
3. integer, intent(IN) :: inod
Process ID

4.17

eigen owner node

Returns owner process ID corresponding to the speciﬁed global loop counter value.
integer function eigen owner node( ictr, nnod, inod )1. integer, intent(IN) :: ictr
Global loop counter
2. integer, intent(IN) :: nnod
Process number
3. integer, intent(IN) :: inod
Process ID

4.18

KMATH EIGEN GEVThis is a generalized eigenvalue computing driver routine that uses EigenExa as the eigenvalue
computing engine. In this driver, eigen_sx is called to compute the eigenpairs via transformation to a pentadiagonal matrix. The constraints on this driver are similar to those on eigen_sx.
Use of this driver routine requires linking KMATH_EIGEN_GEV.o together with EigenExa itself.
subroutine kmath eigen gev( n, a, lda, b, ldb, w, z, ldz )1. integer, intent(IN) :: nMatrix and vector dimensions
2. real(8), intent(INOUT) :: a(lda,*)Matrix A for computing the pencil (A − λB)
The array content is destroyed upon subroutine termination.
3. integer, intent(IN) :: lda
Leading dimension of array a
4. real(8), intent(INOUT) :: b(ldb,*)
Matrix B for computing the pencil (A − λB)
The matrix for transformation to the standard eigenvalue problem is stored
upon subroutine termination.
5. integer, intent(IN) :: ldb
Leading dimension of array b
6. real(8), intent(OUT) :: w(n)Eigenvalues in ascending order
7. real(8), intent(OUT) :: z(ldz,*)
B-orthogonal eigenvectors of the generalized eigenvalue problem
8. integer, intent(IN) :: ldz
Leading dimension of array z

Chapter 5

Other key considerations
5.1

Caution regarding compatibility

As the successor to EigenK, EigenExa has inherited many of its functions. Complete compatibility between the two is not guaranteed, however, since their internal implementations diﬀer
in certain details. These are mainly diﬀerences in function and variable naming rules and incommon domain management methods. For the same reason, simultaneous linking of EigenExa
and EigenK is not recommended.

5.2

Binding with other languages

The method for calling EigenExa from a language other than Fortran90 is highly dependent on
the user’s environment. For further information, refer to “Language bindings” and “Method oflinking to multiple programming languages” in the compiler manual. Information of reference
may also be found in the “Python binding of EigenExa” project [10], which enabled calling fromthe Python language.

5.3

Behavior on error occurrence

During initialization, EigenExa checks that it is being executed under appropriate conditions,but error detection is not performed during execution. In some cases, forced library termination
may occur if a linked subordinate library such as BLAS or LAPACK produces an error.
Information on bug discoveries is essential for improvement of library quality. On discoveryof any bug, please be sure to report it to the developer email address listed on the public websiteof the library.

5.4

Shared library handling in versions 1.x

Shared libraries were not supported in former versions (1.x), because at the time of their development it was not possible to guarantee complete, collision-free resolution of function names
when shared libraries are being used (with certain versions of gcc, abnormal shutdown occurred
without resolution of function names at execution). When version 1.x is to be used as a sharedlibrary, this must therefore be performed solely at the user’s responsibility.
EigenExa versions 2.x and later are shared-library capable, a development achieved withthe technical cooperation of Team Leader Toshiyuki Maeda and other members of the HPC
19

20

CHAPTER 5. OTHER KEY CONSIDERATIONS

Usability Research Team at the RIKEN Advanced Institute for Computational Science. As
noted in the description of the library build, one should select the appropriate make.inc, andwhen executing, always remember to make the appropriate settings for the environment variables(such as LD_LIBRARY_PATH).

Appendix A

Algorithm overview
A.1

Introduction

Appendix A provides an overview of the eigenvalue computation algorithms used in EigenExa,
with the main focus on ontlines of algorithms that two driver routines (eigen s and eigen sx)
use and diﬀerences between them. Both routines are designed to meet the underlying EigenExa
objective of computing all eigenvalues and eigenvectors of real symmetric dense matrices. For
general details about eigenvalue computation algorithms for dense matrices, refer to sources
such as [11, 12, 13, 14, 15].

A.2

Various approaches and related projects
real symmetric (A)ScaLAPACK, EigenExa（（eigen_s））

tridiagonal (T)eigenpairsof Teigenpairsof A

DPLASMA, ELPA
eigenpairsof Bbanded (B)EigenExa（（eigen_sx））

Figure A.1: Various approaches to eigenvalue computation for real symmetric dense matrices.
Let us begin with a brief introduction to the basic aspects of the eigenvalue computing
procedures that is usually applied to real symmetric dense matrices. Textbooks on generalmatrix computation describe an approach based on tridiagonalization of the input matrix (thegreen path in Fig. A.1), which is used in ScaLAPACK [16] (and LAPACK). In the ﬁrst step of
this approach (tridiagonalization), however, the performance is limited by memory bandwidth,
and therefore is expected to be not suﬃciently high on recent computer systems.
This problem led to two development projects, ELPA [17] and DPLASMA [18], which employ an approach based on two-stage tridiagonalization via a banded matrix (the blue path in
Fig. A.1). In this tridiagonalization, the dominant cost arises during the ﬁrst stage, in transformation from dense to band. The byte/ﬂop ratio required in this transformation is smaller
21

22

APPENDIX A. ALGORITHM OVERVIEW

than that in direct tridiagonalization, which means the improvement of eﬀective performance.
But the eigenvector back-transformation process also requires two stages (basically doubling
the cost). Since high-performance implementation in the ﬁrst stage of the back-transformation(from T to B) is currently diﬃcult, its cost becomes extremely large in obtaining a large numberof eigenvectors.
These situations led to the development and provision of two routines for EigenExa. One(eigen s) applies an approach based on the conventional (one-stage) tridiagonalization (the green
path in Fig. A.1). The other (eigen sx) applies an approach in which the eigenvalues andeigenvectors of a banded matrix are computed directly (the red path in Fig. A.1). The following
sections describe these two approaches in a little more detail.

A.3

eigen sAs noted above, the eigen s routine in EigenExa applies an approach based on the conventional(one-stage) tridiagonalization, which is used in ScaLAPACK and other libraries. More specifically, it obtains solutions to the eigenvalue problem Axi = λi xi (i = 1, . . . , N ) through thefollowing three steps:
1. Tridiagonalization of the input matrix by Householder transformations: Q⊤ AQ → T
2. Computation of the eigenvalues and eigenvectors of a tridiagonal matrix by the divideand-conquer method: T y i = λi y i
3. Back transformation of the eigenvectors: Qy i → xi
In step 1, the Householder transformations act from both sides⊤⊤
HN−2 · · · H1 AH1 · · · HN −2 → T,

Hi = I − ui βi u⊤
i

(A.1)

with each column (row) of the input matrix transformed in turn to a tridiagonal matrix(Fig. A.2(a)). Here, we chose the position of the variable beta in the equations for its correspondence with the description further below. The computation of the transform by eachHouseholder transformation is usually performed by using the symmetry of A, as(I − uβu⊤ )⊤ A(I − uβu⊤ ) = A − uv ⊤ − vu⊤ ,1
v = (w − uβ ⊤ (w⊤ u))β,2

w = Au. (A.2)

Furthermore, the Dongarra’s method makes it possible to apply a number of transformations
to the matrix in the form of matrix-matrix multiplications at a time:(I − uK βK uK ⊤ )⊤ · · · (I − u1 β1 u1 ⊤ )⊤ A(I − u1 β1 u1 ⊤ ) · · · (I − uK βK uK ⊤ ) = A − U V ⊤ − V U ⊤ .
(A.3)
However, matrix-vector multiplications, whose performance is limited by the memory bandwidth, still remain (for obtaining the matrix V ) and is therefore a major bottleneck for highperformance computing.
In the second step, the divide-and-conquer method proposed by Cuppen [19] is applied tocompute the eigenvalues and eigenvectors of the tridiagonal matrix. As shown in Fig. A.2(b), a
tridiagonal matrix can be decomposed into a block diagonal matrix and a rank one perturbation.
The underlying idea of the method is computing the eigenvalue decomposition of the tridiagonal
matrix eﬃciently by using the eigenvalue decomposition of the block diagonal matrix (andrecursively apply this idea to the block diagonal matrix).
In the third step, back transformation of the eigenvectors is performed by applying theHouseholder transformations obtained in the ﬁrst step to the eigenvectors of the tridiagonal

A.4. EIGEN SX

23

matrix in the reverse order. Since a number of Householder transformations can be aggregatedinto a convenient form (compact-WY representation):⊤⊤
H1 · · · HK = (I − u1 β1 u⊤1 ) · · · (I − uK βK uK ) → I − U SU ,U = [u1 · · · uK ](A.4)

at low cost (only for computing a small matrix S), the back transformation is usually computed
via matrix-matrix multiplications (Level-3 BLAS):⊤
H1 · · · HN −2 Y = (I − U1 S1 U1⊤ ) · · · (I − UM SM UM)Y → X,

(A.5)

where
Y = [y 1 · · · y N ],X = [x1 · · · xN ].

(A.6)

For this reason, this step is expected to achieve high performance.

=

・・・(a) Tridiagonalization(b)The
method+divide-and-conquer

Figure A.2: Schematic of eigenvalue computation by eigen s.
In eigen s, the ﬁrst and third steps are newly implemented from scratch with appropriate
thread parallelization, whereas the second step is almost ported from the ScaLAPACK code.

A.4

eigen sx

The other driver routine provided in EigenExa, namely eigen sx, is based on an approach that
employs direct computation of the eigenvalues and eigenvectors of a banded matrix. At present,
for the reason mentioned in the last section, a pentadiagonal matrix is selected as the bandedmatrix. More speciﬁcally, the eigenvalue problem is solved in the following three steps.
1. Pentadiagonalization of the input matrix by block version of Householder transformations:Q̃⊤ AQ̃ → B
2. Computation of the eigenvalues and eigenvectors of a pentadiagonal matrix by the divideand-conquer method: By i = λi y i
3. Back transformation of the eigenvectors: Q̃y i → xi
In the ﬁrst step, the block version of Householder transformations are applied to the input
matrix from both sides:⊤
H̃N/2−1· · · H̃1⊤ AH̃1 · · · H̃N/2−1 → P,

H̃i = I − ũi β̃i ũ⊤
i

(A.7)

to transform every two columns (two rows) of the input matrix into a pentadiagonal matrix(Fig. A.3(a)), where
(
)
(i)
(i)
β11 β12(i) (i)ũi = [u1 u1 ], β̃i =
(A.8)
(i)
(i) .
β21 β22
There is no diﬀerence excepting the form of H̃ between Eqs. (A.1) and (A.7), so that the procedure of the pentadiagonalization is the same as the tridiagonalization; the Dongarra’s method
can similarly be applied. The performance bottleneck thus resides in the part of computing Aũ.

24

APPENDIX A. ALGORITHM OVERVIEW

In the second step, as shown in Fig. A.3(b), the pentadiagonal matrix is decomposed into
a block diagonal matrix and a rank two perturbation. By treating the rank two perturbation
as two rank one perturbations, we apply the principal of the divide-and-conquer method fora tridiagonal matrix twice and compute the eigenvalues and eigenvectors of the pentadiagonal
matrix [20].
In the third step, the block version of Householder transformations obtained in the ﬁrst
step are applied to the eigenvectors of the pentadiagonal matrix in the reverse order, whichis essentially the same as in the case of tridiagonalization. The block version of Householdertransformations can also be aggregated in a form with matrices as in Eq. A.4, and matrix-matrix
multiplication can therefore be used in this step, which promises that this step easily achieves
high performance.

=

・・・(a) Pentadiagonalization(b)The
method+divide-and-conquer

Figure A.3: Schematic of eigenvalue computation by eigen sx.
In eigen sx, as in eigen s, steps 1 and 3 are newly implemented with appropriate threadparallelization, whereas step 2 is a simple extended implementation of the ScaLAPACK code
for a pentadiagonal matrix.

A.5

Diﬀerences between eigen s and eigen sx

As described in A.3 and A.4, eigen s and eigen sx comprise three similar steps and are nearly
the same in computation procedures. Particularly in the step of the back transformation of theeigenvectors, there is no essential diﬀerence between them. In this section, we mention the main
diﬀerences in the ﬁrst and second steps between them.
Tri/Pentadiagonalization
The essential diﬀerence between eigen s and eigen sx in this step is that eigen s processes asingle vector, whereas eigen sx processes two vectors together, e.g.
w = Au (in eigen s)→[w1 w2 ] = A[u1 u2 ] (in eigen sx).

(A.9)

In the overall step, the total number of ﬂoating-point operations is a about the same (at least for
the highest term) for the two routines; amount per operation required in eigen sx in about twicethat in eigen s, but the number of operations in eigen sx is about half that in eigen s because
the former deals with two columns at an operation. For the similar reasons, the amount of thedata transfered among distributed processes is almost the same for the two.
The ﬁrst diﬀerence that becomes evident between the two is in the eﬀective performanceof the ﬂoating-point operations (in particular, in matrix-vector multiplications). The data of
matrix A can be reused when computing A[u1 u2 ], whereas it cannot when computing Au.
This means that the required byte/ﬂop ratio in the former is lower than that in the latter. As aresult, the eﬀect of limitation by memory bandwidth is smaller in the former than in the latter(theoretically by about half), which indicates the increasing eﬀective performance in eigen sx.
The second diﬀerence that becomes evident is in the communication latency, arising from thediﬀerence in communication frequencies. Although the data amount per communication is larger

A.6. CONCLUSION

25

in eigen sx than in eigen s, the frequency of communications is lower (by about half). The featureof lower communication frequency (Communication-Avoidance) is a very substantial diﬀerence,
especially in cases of massively parallel computing, due to the fact that communication latency
has become a major problem on recent systems.
In short, eigen sx exhibits clear advantages over eigen s in terms of both the performance of
ﬂoating-point operations and the latency cost of communication. In cases where the problem
size relative to the number of processes is suﬃciently large (with the time for ﬂoating-point
operations thus dominant), the advantage in eﬀective performance is signiﬁcant. On the other
hand, in cases where the number of processes is large (with communication time thus dominant),
the advantage in the latency cost is therefore signiﬁcant. In total, eigen sx is expected to achieve
higher performance than eigen s.
The divide-and-conquer method
It is clear that eigen sx requires more cost (both for ﬂoating-point operations and communications) than eigen s because the former deals with a rank two perturbation whereas the latter
deals with a rank one perturbation. In eigen sx, a rank two perturbation is dealt with two rank
one perturbations. The computational cost for the ﬁrst rank one perturbation can be reduced
by exploiting the structure of the matrix (i.e. block diagonal), however such beneﬁt done notexist in the computation for the second rank one perturbation; the latter cost is about twicethat of the former cost. The resulting cost required in eigen sx increases by a factor of about
three.
In the divide-and-conquer method, one can reduce the cost substantially by the techniqueknown as “deﬂation”. The number of opportunities in which deﬂation can be applied varies
with the problem. In addition, it is diﬀerent even for the same input matrix between routines
via tridiagonalization and pentadiagonalization. Therefore, it is not easy to derive a theoreticalestimation of the diﬀerence between the costs of the two routines.

A.6

Conclusion

In this chapter, we gave an overview of the algorithms employed in the two routines, eigen s andeigen sx, provided in EigenExa and explained their main diﬀerences. Increasing the bandwidthof the banded matrix generally involves the trade-oﬀ; it is advantageous in the ﬁrst step (thetransformation step), but disadvantageous in the second step (the divide-and-conquer method).
Taking this trade-oﬀ into account, we deem that the pentadiagonal matrix is appropriate on
today’s systems. With increasing performance of future systems and improved implementationof the divide-and-conquer method (our ScaLAPACK-based implementation seems to be rarely
suitable for current systems), banded matrix with larger band width (e.g. heptadiagonal) may
prove promising. By contrast, the use of conventional tridiagonalization (eigen s) might in somecircumstances be the best choice. We hope that an understanding of the information in thisappendix will help users to select the routine that is most appropriate to their application.

26

APPENDIX A. ALGORITHM OVERVIEW

Acknowledgements
We would like to express our deepest gratitude to all the members of the EigenExa developmentgroup for their sincere contribution to developing EigenExa. We also grateful the support fromthe RIKEN Advanced Institute for Computational Science (AICS). In addition, we appreciate
the users of EigenExa for their bug reports and feedback that is beneﬁcial to further quality
improvement of EigenExa. Without people and support above mentioned, EigenExa could havenot been released.
The EigenExa project has been supported by the following funds.
• JST CREST, “Development of System Software Technologies for post-Peta Scale High
Performance Computing” (FY2011–FY2015).
• MEXT KAKENHI, Grant Numbers 21300013 (FY2012), 23240005 (FY2011–FY2013), and
15H02709 (FY2015–FY2017).
The EigenExa project has used computational resources of the K computer provided by the
RIKEN AICS through the folowing projects.
• HPCI System Research project, project IDs hp140069 (FY2014) and hp140069 (FY2012–
FY2013).
• Computational resources for enhancement, Project ID ra000005 (FY2013–).
The EigenK project, the predecessor of the EigenExa project, was supported by the following
fund.
• JST CREST, “High Performance Computing for Multi-Scale and Multi-Physics Phenomena” (FY2006–FY2012).

27

28

APPENDIX A. ALGORITHM OVERVIEW

References[1] S. Yamada, T. Imamura, T. Kano and M. Machida, “High-Performance Computing for
Exact Numerical Approaches to Quantum Many-Body Problems on the Earth Simulator”,
Proceedings of the 2006 ACM/IEEE Conference on Supercomputing (SC ’06), November
2006. Tampa USA.
http://doi.acm.org/10.1145/1188455.1188504[2] T. Imamura, S. Yamada and M. Machida, “Development of a High Performance Eigensolver
on the Peta-Scale Next Generation Supercomputer System”, Progress in Nuclear Scienceand Technology, the Atomic Energy Society of Japan, Vol. 2, pp.643–650 (2011) .
[3] T. Imamura, S. Yamada and M. Machida, “Eigen-K: high performance eigenvalue solver forsymmetric matrices developed for K computer”, 7th International Workshop on. Parallel
Matrix Algorithms and Applications (PMAA2012), June 2012, London UK.
[4] What is K? — RIKEN Advanced Institute for Computational Science,
http://www.aics.riken.jp/en/k-computer/about/[5] K computer — Fujitsu Global,
http://www.fujitsu.com/global/about/businesspolicy/tech/k/[6] T. Imamura and Y. Yamamoto, “CREST: Dense Eigen-Engine Groups”, International
Workshop on Eigenvalue Problems: Algorithms; Software and Applications, in PetascaleComputing (EPASA 2014), Tsukuba, March 7–9, 2014 (poster).
http://www.aics.riken.jp/labs/lpnctrt/EPASA2014_dense_poster_ImamuraT_only.
pdf[7] T. Imamura, “The EigenExa Library – High Performance & Scalable Direct Eigensolver forLarge-Scale Computational Science”, HPC in Asia, Leipzig, Germany, June 22–26, 2014.
[8] T. Imamura, Y. Hirota, T. Fukaya, S. Yamada and M. Machida, “EigenExa: high performance dense eigensolver, present and future”, 8th International Workshop on Parallel
Matrix Algorithms and Applications (PMAA14), Lugano, Switzerland, July 2–4, 2014.
[9] T. Fukaya and T. Imamura, “Performance evaluation of the EigenExa eigensolver on the
Oakleaf-FX supercomputing system”, Annual Meeting on Advanced Computing System
and Infrastructure (ACSI) 2015, Tsukuba, January 26–28, 2015.
[10] Python binding of EigenExa, HPC Usability Research Team, RIKEN AICS, http://www.
hpcu.aics.riken.jp/[11] B. Parlett, “The Symmetric Eigenvalue Problem”, SIAM (1987).
[12] J. Demmel, “Applied Numerical Linear Algebra”, SIAM (1997).
[13] L. Trefethen and D. Bau, III, “Numerical Liner Algebra”, SIAM (1997).
29

30

REFERENCES[14] W. Press, S. Teukolsky, W. Vetterling and B. Flannery, “Numerical Recipes: The Art of
Scientiﬁc Computing”, 3rd ed., Cambridge University Press (2007).
[15] G. Golub and C. Van Loan, “Matrix Computations”, 4th ed., Johns Hopkins University
Press (2012).
[16] ScaLAPACK, http://www.netlib.org/scalapack/[17] ELPA, http://elpa.rzg.mpg.de/[18] DPLASMA, http://icl.cs.utk.edu/dplasma/[19] J. Cuppen, “A divide and conquer method for the symmetric tridiagonal eigenproblem”,
Numer. Math, Vol.36, pp.177–195 (1981).
[20] P. Arbenz, “Divide and conquer algorithms for the bandsymmetric eigenvalue problem”,
Parallel Computing, Vol.18, No.10, pp.1105–1128 (1992).