Part I
Architecture and Algorithms ofHigh Performance Computing

High-speed ultra-low-power simulation methodology
for material science
Mary Inaba1,2, Hiroshi Imai1,2, Reiji Suda1,2, Kei Hiraki1,3
The University of Tokyo1
Co-Investigator(Kenkyu-Buntansha)2
Co-Investigator(Renkei Kenkyusha)3

1. Outline	  2011 - 2014 Research Project
We performed studies on High-speed ultra-low-power simulation methodology formaterial science. The main goal of the study is to establish technologies forhigh-speed computation, high-speed network and high-speed memory systems toaccelerate simulations in material science area. Our focus is the use of FPGA (field
programmable gate array) for accelerating computation, network operation andmemory access acceleration.
During 2010 – 2014 periods, we have constructed two computation accelerator, two
network accelerator and several other FPGA board to perform network and memory
emulation.

As results of these studies, we have established basic architecture for

the FPGA based platform for simulation in wide application area including materialscience.
Another important contribution of this project is establishing a new educational
unit, “Computing Science Educational Course” as a joint activity between thedepartment of computer science, graduate school of information science andtechnology and the department of physics, graduate school of science at the
University of Tokyo. This educational course is a result of collaboration in thisproject.

2.

Research topics in 2010 – building basic structure for the study.

2010 was the starting year of this project. We compiled the R&D team, then we

1

performed

study

on

measurementsofprocessor

architecture

and

power

consumption, constructing FPGA based computation accelerator, high-efficiency
data-transfer on long distance high-speed internet, high-speed high-productivity
language HPC Ruby and fast GPGPU-host data-transfer mechanisms.
Figure 1 shows results of GPGPU-host data-transfer optimization, where parallel
data-transfer is scheduled to minimized waiting time of the PCIE bus connecting
GPGPUs and a host. Evaluation is performed using 8192×256 matrix multiply
with Nvidia GPGPU.

Fig 1. Optimization of GPGPU-host data-transfer
3. Research topics in 2011 – FPGA based accelerator and HPC Ruby
We developed a computation accelerator that use very large FPGA that have both a
PCIE interface and a 40Gbps Ethernet interface (Figure 2). Very high-speed
Ethernet connection (40Gbps) gives very short connection latency between
computation accelerators (less than 1 μsecond) compared to a large latency using
PCIE (about 10 μsecond including host drivers). Short latency connections are
essential to realize scientific applications that have strong-scaling computation
behavior.
In 2011 we also develop HPC ruby system that consists of static HPC rubycompiler, type inference system and runtime systems. Since Ruby is a dynamic
weakly typed language with very-high productivity, the only shortcoming of Ruby
was its very slow execution speed. The reason why original Ruby is so slow is that
Ruby is a weakly typed dynamic language that require runtime analysis of data
type and variable bindings. Figure 3 shows speed-up ratio of HPC-Ruby vs. Ruby
1.8 on various numerical application programs. This figure shows clear advantage
of static analysis in dynamic language for numerical computations.

2Fig 2. Block diagram of Computation

Fig 3.

Accelerator

Speedup ratio of HPC Ruby fromoriginal Ruby 1.8

4,Research topics in 2012 – Research on deficiency of compilers and processors
We investigate requirements for future HPC processor that aims EXA flops bymeasuring power performance of SPARC64 XIfx, Intel IvyBridge, IBM Power7 and
IBM BlueGene/Q processor. The results show large difference of absoluteperformance and power efficiency. In order to investigate the performancedifference in detail, we compared performance of programs using several different
compilers (Intel C compiler, GCC and Fujitsu C compiler).
Figure 4 shows performance difference in SPEC CPU2006flt performance of
SPARC64 IXfx with FCC and Intel IvyBridge Xeon with Intel compiler. Due to thelimitation of the space we omit results from other benchmarks and othercombination of processors and compilers, these results shows large performancedifference even if they are normalized by CPU clock.
The largest reason of performance difference is integer and memory performanceof processors. As shown in figure 5, even in floating point heavy applications, more
integer and memory instructions than floating point arithmetic instructions are
executed in many applications. As a result, integer and memory performance ismore important than floating point performance especially in simulations inmaterial science where the problems are less regular than commonly used CFD or
FEM execution.
We further investigate the cause of integer and memory performance. We foundthat part of integer/memory performance is caused by the difference in optimizationof compiled codes.

3Figure 4. Performance comparison of SPARC64 IXfs with Fujitsu C compiler andIntel Xeon E5-2687 with Intel C compiler . Benchmark programs are SPEC
CFP2006.

size=S	
 
ua
mg
is

Int

ep

Fp

bt
0%

Figure 5.

20%

40%

60%

80%

100%

Percentage of integer and memory instructions executed in NASPAR

benchmarks with Intel C compiler.
We compared the performance and object codes of (1) GCC 4.5.3, (2) Intel Composer XE
2013 and Fujitsu C compiler V1.2.1.
For example, we compiled the source code as:
Proc_3 (Rec_Pointer *Ptr_Ref_Par){ if (Ptr_Glob != Null)*Ptr_Ref_Par = Ptr_Glob->Ptr_Comp;4

Proc_7 (10, Int_Glob, &Ptr_Glob->variant.var_1.Int_Comp);}
FCC outputs 22 instructions, GCC outputs 15 instructions. The main performance andcode size difference between these two are-

Unnecessary save/restore-

Partial use of available addressing modes of instructions-

Address calculation is performed by combination of general-purpose instructions-

FCC does not use brz branch instruction. Instead it use a combination of cmp and be

These absence of basic optimization cancelled many complicated optimization include in
FCC.
Consequently, INT performance / power consumption becomes worse than expected.

Figure 6.

Yearly progress in INT/memory processing speed/W
Red line shows that of SPARC64 IXfx (FX10)

Figure 6 shows yearly progress in INT/memory processing speed /W by executing
Dhrystone benchmarks for various processors. This figure shows about 8 times
MIPS/W improvement is necessary to catch up latest efficient processors, and morethan 16 times is necessary to construct Exa Flops system only with general purpose
processors.
５．Research topics in 2013 – FPGA base network accelerator prototype
In 2015, we developed 2nd

FPGA based network accelerator that use large-scale

FPGA with ARM processor cores (Figure 7). Interconnections between FPGA boards

5are realized by 24-port SATA 6Gbps connections that can be dynamicallyre-configured. This FPGA based network emulator have been used to evaluatevarious design parameters of an interconnection network in a parallel system andnetwork on chip (NoC) for many-core processor systems.
Figure 8 shows an example of resource evaluation for various network topologies.
As shown in Figure 8, mesh, ring, torus and shuffle-exchange network have almost
the same amount of FPGA hardware resource because major part of the resource isused to make router buffers.

Figure 7.

FPGA based network accelerator and its block diagram

Figure 8 FPGA hardware budgets of different network topologies
6. Research topics in 2014 -- Optimization of hardware resource for networks
During 2014, we conducted studies on (1) Network performance optimization using
Genetic Algorithm, (2) FPGA based network accelerator/emulator (continued from2012,2013), (3) reduction of power consumption of networks by dynamic routing
function,(4) FPGA based computation accelerator for machine learning by neural6

networks and (5) Software DSM to access PGAS with LRC protocol.

Fig 9. Optimized to random traffic

Fig 10. Optimized to transpose

Figure 11. Torus network optimized for various traffic
Due to space limitation, we showed results of Network performance optimization
using Genetic Algorithm.
Almost all the interconnection network, especially network on chip(NoC) has been
designed manually without optimization of the width of network links and capacityof routing nodes. As a result, performance and power consumption of networks are
not fit to their major applications.
We adopt genetic algorithm (GA) to optimize heterogeneous networks thatmaximize performance and power efficiency. In this study, all the links can beselected from fast links and slow links where the bandwidth is 2 times different.
Also, all the routers can be selected from large (buffer) routers and small routers.
We applied GA to find the most optimized configuration of the network to transfer acertain defined traffic patterns such as random, transpose, all to all etc. We perform
two optimization with GA, where the one assumes X and Y symmetric property ofthe target traffic,, and the one without this assumption.
The target topology of the network is either 2-D mesh or 2D-torus. Figures 8 and 9is the result for 2-D mesh and Figure 10 is for 3-D torus.	
  These results show that
equally distributed mesh or torus is far worse from optimized configuration. Thismeans, commonly used mesh network or torus network of7

many core processors

are not optimal and gives worse power consumption.
We are now performing more detailed evaluation using real network traffic ofmajor application programs.
Selected Publications 2010 -- 2015[1] Hisanobu Tomari, Mary Inaba, Kei Hiraki, “Compressing Floating-Point
Number Stream for Numerical Applications”, The First International Conference on
Networking and Computing, Hiroshima, 2010.
[2] Seongsoo Moon, Mary Inaba, “Modeling and simulation for EMS”, 3rdInternational

Conferenceon

Wireless

Information

Networks&

Business

Information System (WINBIS'11), pp.94-100, Kathmandu Nepal, 2011.
[3] T. Tanuma, H. Imai and S. Moriyama, “Revisiting Hyperbolic Voronoi Diagramsfrom Theoretical, Applied and Generalized Viewpoints”, 2010 International
Symposium on Voronoi Diagrams in Science and Engineering (ISVD '10), pp.23-32,
Quebec, 2010.
[4] Kamil Rocki, Reiji Suda, “MPI-GPU Monte Carlo Tree Search”, IEEE 2011
International Conference on Information and Computer Applications(ICICA), Dubai,
UAE, 2011.
[5] Kamil Rocki and Reiji Suda, “Massively Parallel Monte Carlo Tree Search”,
Proceedings of VECPAR 2010, 9th International Meeting High Performance
Computing for Computational Science, 2010.
[6] Yasuo Ishii, Mary Inaba, Kei Hiraki, “Access Map Pattern Matching for High
Performance Data Cache Prefetch”, The Journal of Instruction Level Parallelism,
Vol.13, (ISSN 1942-9525), pp.1-24, 2011.
[7] R. Suda, K. Naono, K. Teranishi and J. Cavazos, “Software Automatic Tuning:
Concepts and State-of-the-Art Results”, Chapter 1 of “Software Automatic Tuning:
From Concepts to the State-of-the-Art Results”, Springer, Sep. 2010. pp.3-18.
[8] R. Suda, “A Bayesian Method of Online Automatic Tuning”, Chapter 16 of“Software Automatic Tuning: From Concepts to the State-of-the-Art Results”,
Springer, Sep. 2010.
[9] Junichiro Makino, Hiroshi Daisaka, Toshiyuki Fukushige, Yutaka Sugawara,
Mary Inaba, Kei Hiraki, “The performance of GRAPE-DR for dense matrix
operations”, the International Conference on Computational Science, ICCS 2011,
Nanyang, Singapore, 2011.
[10] Yasuo Ishii, Keisuke Kuroyanagi, Takeo Sawada, Mary Inaba, Kei Hiraki,“Revisiting Local History for Improving Fused Two-Level Branch Predictor”, 2nd

8

JILP Workshop on Computer Architecture Competitions (JWAC-2): Championship
Branch Prediction, San Jose, USA, Jun4, 2011.
[11] Tomohiro Sonobe, Mary Inaba, Ayumu Nagai, “Counter Implication Restart”,
Pragmatics of SAT 2011,

Ann Arbor, USA, Jun.18, 2011[12] Tomohiro Sonobe, Mary Inaba, “Counter Implication Restart for Parallel SAT
Solvers”, Learning and Intelligent Optimization Conference (LION6), Paris, France,
2012[13] Yasuo Ishii, Kouhei Hosokawa, Mary Inaba, and Kei Hiraki,“High Performance
MemoryAccessSchedulingusing

Compute-Phase

Prediction

andWriteback-Refresh Overlap”, 3rd JILP Workshop on Computer Architecture
Competitions (JWAC-3): Memory Scheduling Championship (MSC), 2012.
[14] Tomohiro Sonobe and Mary Inaba,“Division and Alternation of Decision
Variables”, Pragmatics of SAT 2012, 2012.
[15] Yasuo Ishii, Mary Inaba, Kei Hiraki,“Unified Memory Optimizing Architecture:
Memory Subsystem Control with a Unified Pedictor”, 26th International
Conference on Supercomputing, 2012.
[16] Tomohiro Sonobe and Mary Inaba,“Portfolio with Block Branching for Parallel
SAT Solvers”, Learning and Intelligent OptimizatioN conference (LION7), 2012.
[17] Toru Motoya and Reiji Suda, “Conjugate Gradient Methods Relieved for Inner
Product Communication Latencies”, International workshop on HPC, Krylov
Subspace method and its applications, 2012.
[18]

Reiji

Suda,“Toward

TunableMulti-Scheme

Parallelization”,

SIAM

Computational Science and Engineering, 2013.
[19] V. Suppakitpaisarn, M. Edahiro, and H. Imai,“Fastest Multi-Scalar
Multiplication Based on Optimal Double-Base Chains”, World Congress on Internet
Security (WorldCIS-2012), 2012.
[20] N. Fu, A. Hashikura and H. Imai,“Geometrical treatment of periodic graphswith coordinate system using axis-fiber and an application to a motion planning”,
The 9th Int. Symposium on Voronoi Diagrams in Science and Engineering, 2012.
[21] C.Chuengsatiansup, H.Imai, and V.Suppakitpaisarn,“Evaluating Optimized
Computation in Double-Base Chain for Efficient Elliptic Curve Cryptography”, 15th
Japan-Korea Joint Workshop on Algorithms and Computation, 2012.
[22] V. Suppakitpaisarn, M. Edahiro, and H. Imai,“Fast Elliptic Curve
Cryptography Using Optimal Double-Base Chains”, the International Conferenceon Informatics & Application (ICIA2012), 227, pp190-204, 2012.
[23] V. Suppakitpaisarn, M. Edahiro, and H. Imai,“Optimal Elliptic Curve Scalar

9Multiplication Using Double-Base Chains”, International Journal of Digital
Information and Wireless Communications (IJDIWC), Vol. 2, pp115-134, 2012.
[24] T. Satoh, F. Le Gall, and H. Imai,“Quantum network coding for quantum
repeaters”, Phys. Rev. A 86, 032331, 2012.
[25] Norie Fu, Akihiro Hashikura, and Hiroshi Imai,“Proximity and Motion
Planning on 1-Rigid Planar Periodic Graphs”, Transactions on Computational
Science, Vol.XX (LNCS Vol.8110), pp.72–91, 2013.
[26] Reiji Suda, Luo Cheng, Takahiro Katagiri,“A Mathematical Method for Online
Autotuning of Power andEnergy Consumption with Corrected Temperature
Effects”Procedia Computer Science, Vol 18, pp 1302-1311, ISSN 1877-0509, 2013.
[27] Kei Hiraki, Mary Inaba, Akira Kato, “Efficient utilization of 40/100 Gbps
long-distance network by a single TCP stream”, TERENA NETWORKINGCONFERENCE (TNC2013), 2013.
[28] Tomoaki Ishiyama, Junichiro Makino, Simon Portegies Zwart, Derek Groen,
Keigo Nitadori, Steven Rieder, Cees de Laat, Stephen McMillan, Kei Hiraki, Stefan
Harfst,“The Cosmogrid Simulation: Statistical Properties of Small Dark Matter
Halos” The Astrophysical Journal, DOI: 10.1088/0004-637X/767/2/146.
[29] J.-F. Baffier, V. Suppakitpaisarn, H. Hiraishi and H. Imai, “ParametricMultiroute Flow and Its Application to Robust Network with k Edge Failures”,
Lecture Notes in Computer Science, Vol.8596, pp.26-37 (2014).
[30] H. Hiraishi and H. Imai, “BDD Operations for Quantum Graph States”,
Lecture Notes in Computer Science, Vol.8507, pp.216-229 (2014).
[31] Tomohiro Sonobe, Shuya Kondoh and Mary Inaba, “Community Branching for
Parallel Portfolio SAT Solvers”, The International Conferences on Theory andApplications of Satisfiability Testing (SAT2014) pp. 188-196, 2014.
[32] Hongzhi Chen, Reiji Suda, “A Noise-Reducing Collective Communication
Algorithms”, 2015 Conference on Advanced Topics and Auto Tuning in
High-Performance Scientific Computing (ATAT HPSC 2015), National Taiwan
University, Feb. 27-28 (27), 2015.
[33] Takahiro Naruko, Kei Hiraki, “Reducing Energy Consumption of NoC by
Energy Efficient Router Bypassing (EERB),” Transaction on ACS, Information
Processing Society of Japan, ACS）,8(1),12-21, 2015.
[34] Jun Ohno, Kei Hiraki, “Accelerating cache coherence mechanism with
speculation,” poster

ICS 2014: 171, 2014[35] Tomohiro Sonobe, Shuya Kondoh, Mary Inaba, “Community Branching for
Parallel Portfolio SAT Solvers,”

SAT 2014: 188-196, 201410

Numerical Computation Algorithms for Large-scale Parallel Environment
Daisuke Takahashi1, Toshiyuki Imamura2, Hiroto Tadano1
University of Tsukuba1, RIKEN, Advanced Institute for Computational Science2
Research was conducted on the use of large-scale parallel environments, whose
performance exceeded 10 petaflops, for numerical computation algorithms. These
included the fast Fourier transform (FFT), quadruple- and octuple-precision basic linear
algebra subprograms (BLAS) on graphics processing units (GPUs), triple- andquadruple-precision floating point operations on GPUs, sparse matrix-vector
multiplication on GPUs, and a real symmetric eigenvalue solver operating in a CPU +GPU environment. An evaluation was also carried out of the performance of a newly
developed eigenvalue solver and an existing solver in a GPU environment, ahigh-performance eigenvalue solver employing third-generation NVIDIA GPUs, solution
acceleration of simultaneous linear algebraic equations using the Block Krylov algorithm,
accuracy improvement of Block Krylov partial space iterative methods, andperformance-tuning methods.

1. Introduction
As of March 2015, five supercomputer systems (one of which incorporate GPUs) haveappeared with performances exceeding 10 petaflops. Future mainstream supercomputers are
expected to incorporate multi-core CPUs as well as accelerators such as GPUs as computation
nodes, of which there will be several thousand or several tens of thousands.
However, supercomputers whose configurations incorporate multi-core CPUs (and multipleGPUs) are beginning to experience difficulty in gaining high execution efficiency as aconsequence of the increased processor core count and deficient memory bandwidth per
operational performance unit. Therefore, the parallel numerical computing algorithms andperformance-tuning methods proposed until now will be insufficient to meet futurecomputational science challenges.
The following research topics have therefore been selected in the present study, andperformance-tuning methods were also investigated:- Parallel three-dimensional Fast Fourier transform (FFT) algorithm for next-generation
supercomputers- Performance evaluation of quadruple- and octuple-precision BLAS implementation on
GPUs- Acceleration of sparse matrix-vector multiplication by automatic storage format selection
with GPU- Performance evaluation of parallel two-dimensional FFT using Intel AVX instructions- Study of triple- and quadruple-precision floating point operations with GPU- Optimization and performance tuning technique for a numerical software on Multi-{core,socket, GPU} computing environments- Development of efficient Block Krylov subspace methods for linear systems with multiple
right-hand sides and their applications
These parallel numerical computing algorithms and performance-tuning methods are11expected to allow computation on a scale previously considered impossible when they areapplied to actual programs used in computational material science.

2. Parallel three-dimensional FFT algorithm for next-generation supercomputers
A typical decomposition for performing a parallel three-dimensional FFT is slab-wise. In thiscase, a three-dimensional array x(Nx, Ny, Nz) is distributed along the third dimension Nz.
However, Nz must be greater than or equal to the number of MPI processes. This becomes anissue with very large node counts for clusters.
To solve this problem, a scalable framework for three-dimensional FFTs on Blue Gene/L
supercomputer and a three-dimensional FFT on the six-dimensional network torus QCDOC
parallel supercomputer have been proposed. The parallel three-dimensional FFT algorithms use
the three-dimensional decomposition for the three-dimensional FFT. These schemes require
three all-to-all communications.
We have proposed a two-dimensional decomposition for the parallel three-dimensional FFT.
This requires only two all-to-all communications.
In addition to reducing the communication time through two-dimensional decomposition,
overlapping operations and communications further improved the performance compared toprevious implementations.
The parallel three-dimensional FFT performance for the T2K-Tsukaba is shown in Fig. 1.
Performance improvement through overlapping of operations and communications is evident.

Fig. 1. Parallel three-dimensional FFT performance for T2K-Tsukuba(64 nodes, 256 MPI processes × 4 threads)

2. Performance evaluation of quadruple- and octuple-precision BLAS implementation on GPUs
BLAS functions supporting quadruple- and octuple-precision operations were implemented
on GPUs and evaluated. For quadruple- and octuple-precision operations, 2 double-type
double-precision numbers were linked to represent quadruple-precision numbers indouble-double (DD) type quadruple precision operations, and 4 double-type double-precision
numbers were linked to express octuple-precision numbers in quad-double (QD) type octuple
precision operations.
Performance evaluation with an NVIDIA Tesla C2050 revealed that quadruple-precision
AXPY was accelerated approximately 9.5 times and octuple-precision AXPY was accelerated

12

approximately 19 times compared to an Intel Core i7 920 for the same operations. In addition,
quadruple-precision general matrix multiplication (GEMM) was accelerated approximately 29
times compared to the CPU, and octuple-precision GEMM was accelerated approximately 24
times. The quadruple-precision AXPY for the Tesla C2050 was only 2.1 times the
double-precision operations, and the increased computing times for double-precision operations
with the GEMV and GEMM were also greatly reduced in comparison to the CPU case.
When PCI-Express (PCIe) data transfer times were considered, double-precision GEMM
showed a tendency to be rate-limited by the PCIe data transfer rate, but quadruple- andoctuple-precision GEMM essentially eliminated this limit.

3. Acceleration of sparse matrix-vector multiplication by automatic storage format selection
with GPU
Sparse matrices are frequently the target of simultaneous linear algebraic equations inscientific and engineering computations. Solutions are derived in this case through iterative
methods, and thus acceleration of sparse vector-matrix multiplication is desired. Various storage
formats for sparse matrices have been proposed. Depending on the target sparse matrix,
performance differences emerge for sparse matrix-vector multiplication. Since sparse
matrix-vector multiplication involves more memory access than computation operations,
optimum storage format selection is important for enhancing GPU performance.
This research focused on automatic selection of the storage format for sparse matrices inorder to enable acceleration of sparse matrix-vector multiplication. As a preliminary assessment,
sparse matrix-vector multiplication was implemented with four types of storage formats for
performance evaluation. Consequently, performance differences were evaluated based on
variances in non-zero elements and percentages of non-zero elements. Based on the results, analgorithm was implemented for the automatic selection of the optimum storage format, andsolutions of simultaneous linear algebraic equations employing iterative methods were used forevaluation.
According to the performance evaluation results with an NVIDIA Tesla C2050, optimum
storage format selection occurred for many of the sparse matrices employed, and enabled speedsthat were comparable to or higher than that achieved using CUSPARE directly with the CRS
format. In addition, the proposed automatic selection was found to be effective even for very
low numbers of iterations.

4. Performance evaluation of parallel two-dimensional FFT using Intel AVX instructions
Recent processors such as the Intel Xeon SSE, SSE2, SSE3, SSSE3, SSE4, and AVX; theMotorola PowerPC AltiVec; and the Fujitsu SPARC64 Viiifx HPC-ACE often incorporate Short
Vector SIMD commands to accelerate processing of floating point operations.
The data supply performance of recent processors, however, relies on cached memory. Even
if Short Vector SIMD commands are used, memory access optimization needs to be conducted.
The FFT library FFTE was implemented to improve the performance of the FFT kernel portion
by employing AVX commands, and by using a six-step FFT block algorithm, the highperformance of FFTE was maintained even when the data could not fit within the cachedmemory.
The parallel two-dimensional FFT performance is shown in Fig. 2 for the Intel Xeon E5-2670(Sandy Bridge 2.6 GHz, 20MB L3 cache, 2 CPUs, 16 cores, 128 GB DDR3-SDRAM). It can beseen that the performance of FFTE is higher than that of FFTW.

13

Fig. 2. Parallel two-dimensional FFT performance with Intel Xeon E5-2670(Sandy Bridge-EP 2.6 GHz, 20 MB L3 cache, 2 CPUs, 16 cores)

5. Study of triple- and quadruple-precision floating point operations with GPU
In recent years, memory and network bandwidth shortages have become an issue as processor
performance improves. When precision was insufficient in floating point operations withdouble-precision operations, the employment of quadruple-precision operations had beenstudied previously, but cases are believed to exist where triple-precision operations with lessdata-access volume can be valid. In the present study, triple-precision numbers were treated as a
combination of double and single numbers, which were stored separately (D+S type).
Performance was evaluated by implementation of a triple-precision BLAS routine on a Tesla
M2090 GPU. D+S-type operations implemented a method in which some operations were
performed as single-precision operations for the double-double type quadruple-precision
operation (DD-type operation) algorithm. Type conversion between double-precision numbers
and single-precision numbers occurred frequently, however, and D+S-type operations became
higher in cost than DD-type operations. Consequently, a system in which BLAS I/O wasconducted with the D+S type, and computing using DD-type operations was implemented.
During performance evaluation on the Tesla M2090, triple- and quadruple-precision AXPY and
GEMV became memory rate-limited, and their execution times were shown to be respectively 3
and 4 times more strongly proportional to data size than single-precision routines, and 1.5 and 2
times more strongly proportional than double precision routines.
Although GPUs have a high floating-point operation performance, PCI Express and network
bandwidths are far lower than that for global memory, and can frequently become performancebottlenecks in a GPU cluster environment. The memory space on a GPU is limited, and thislimitation can cause communications with the host via PCI Express. In such environments,
cases are believed to exist where triple precision with smaller data sizes than quadruple
precision is valid.

146. Optimization and performance tuning technique for a numerical software on Multi-{core, socket,
GPU} computing environments
6.1. Introduction
In this study, we focused on the potential of GPGPU, which have been widely spread over not
only HPC but general purpose computing, for running a highly advanced material science programcode. The essence of the study is to accelerate an eigenvalue solver on a desktop class computer
system in order to contribute to this research program shortly. We mainly studied two topics in thisprogram. First, we accelerated the performance of CUDA BLAS library by introducing a novel
CUDA GPU algorithm. Second, we developed a new GPU-based eigenvalue solver by taking
advantage of newly developing CUDA-BLAS functions. This report briefly summarizes theirresults.

Fig. 3. High-end GPU card
(nVIDIA Tesla K20c, peak performance 1.17TFLOPS/Double Precision)
6.2. New Atomic-operation-based SYMV algorithm on a CUDA GPU
In literatures [8,9,18], we studied a new SYMV (symmetric matrix vector multiplication routine)
implementation by taking advantage of utilizing the atomic operations. In the study, ahigh-performance SYMV kernel was implemented on several GPU generation cores using an
atomic-operation based algorithm. The algorithm is effective for the memory bandwidth and reduced
memory usage. In fact, the latest implementation demonstrates sustained double-precision andsingle-precision performances of approximately 67 GFLOPS and 130 GFLOPS, respectively, on a
Tesla K20c GPU card [18]. This is also nearly peak performance of the memory bandwidth, and thealgorithm outperforms major CUDA BLAS kernels; CUBLAS, KBLAS, MAGMABLAS, and
CULA-BLAS. By this novel algorithm, we can expect extreme performance improvement on a GPU
computer system since eigenvalue computations in most of material simulation code bounds bymemory bottleneck. It has a significant impact to the simulation code on which the SYMV kernel isplugged into.

15

Fig. 4. Performance evaluation of the DSYMV kernels (for lower triangular format) in the latest
CUDA-BLAS libraries (ASPEN.K2 (‘+’ marker) is our development).
6.3. GPU-based eigenvalue solver
In literatures [14,3,19], we studied a GPU-based eigenvalue solver by taking advantage of newly
developed SYMV kernel presented in the previous topic. We developed Eigen-G, which is a
GPU-based eigenvalue solver for real-symmetric matrices. We confirmed that Eigen-G outperforms
state-of-the-art GPU-based eigensolvers such as magma dsyevd and magma dsyevd 2stage
implemented in the MAGMA version 1.4.0, when we select appropriate CUDA kernels called in
the solver. In fact, we apply the best-performing CUDA BLAS libraries and the GPU-CPU hybrid
DGEMM technologies (see Fig. 5) into the Eigen-G, the Eigen-G solver yields an even better
performance improvement. In the latest report of Eigen-G [19] reveals good performance
improvement (2x to 3x) compared with the CPU-based eigenvalue solver. We believe that thisperformance improvement must give a big impact to the material science simulation code.

Fig. 5. Two Typical GPU offload model studied in Eigen-G16

7. Development of efficient Block Krylov subspace methods for linear systems with multiple
right-hand sides and their applications
Linear systems with multiple right-hand sides:

appear in many scientific applications such as lattice quantum chromodynamics (QCD)
calculation, an eigensolver based on contour integral, and so on. As numerical solvers forsolving these linear systems, Block Krylov subspace methods have been proposed. However,
the relative residual norm of Block Krylov subspace methods may not converge due to the
influence of numerical instability when the number L of right-hand sides is large. In this project,
we improve the convergence behavior and accuracy of approximate solutions of Block Krylov
subspace methods.
7.1. Improvement of the convergence behavior of the Block BiCGSTAB method
As one of Block Krylov subspace methods, the Block BiCGSTAB method has been proposed.
The number of iterations of the Block BiCGSTAB method may smaller than that of the
BiCGSTAB method. However, if column vectors of n×L matrices of the Block BiCGSTAB
method are numerically linearly dependent, the method becomes unstable. We improve theconvergence behavior of the Block BiCGSTAB method by performing orthonormalization ofcolumn vectors.
Fig. 6 shows the relative residual history of the Block BiCGSTAB method and the Block
BiCGSTAB method with orthonormalization. Test problem is a linear system derived fromlattice quantum chromodynamics calculation. The size n of the coefficient matrix is 1,572,864,
and the number L of the right-hand sides is 8. From Fig. 6, the relative residual norm of theBlock BiCGSTAB method diverged. On the other hand, that of the Block BiCGSTAB methodwith orthonormalization converged.

Fig. 6. Relative residual history of the Block BiCGSTAB method and the Block BiCGSTAB
method with orthonormalization. ■: Block BiCGSTAB method，■: Block BiCGSTAB methodwith orthonormalization.
7.2. Accuracy improvement of approximate solutions of the Block BiCR method
The Block BiCR method is a natural extension of the BiCR method for linear system withsingle right-hand side. This method often shows smooth convergence behavior compared with
the Block BiCG method. However, approximate solutions generated by the Block BiCR method17may deteriorate due to an error matrix arise from the matrix multiplication. We modify thealgorithm of the Block BiCR method to improve accuracy of approximate solutions.
The test problem is Si5H12 from The University of Florida Sparse Matrix Collection. The
size n of the matrix is 19,896, and the number L of the right-hand sides is 32. Fig. 7 shows thetrue relative residual history of the Block BiCG method, the Block BiCR method, and themodified Block BiCR method. In this numerical experiment, residual orthonormalization isapplied to improve numerical stability. The true relative residual norm is used for evaluating
accuracy of the obtained approximate solution. If the true relative residual norm is sufficiently
small, accuracy of the obtained approximate solution is good. From Fig. 7, the true relative
residual norm of the Block BiCR method stagnates around 3.0×10-11. On the other hand, that ofmodified Block BiCR method decreases to 4.3×10-13.

Fig. 7. True relative residual history of the Block Krylov subspace methods. ■: Block BiCG
method，■: Block BiCR method，■: Modified Block BiCR method.
7.3. Speed-up of the eigensolver based on contour integral by using Block Krylov subspacemethods
The eigensolver based on contour integral (SS method: Sakurai-Sugiura method) has been
proposed. The SS method can compute eigenvalues in a specified region on the complex plane
and their corresponding eigenvectors. In this method, linear systems with multiple right-handsides need to be solved. The solution of these linear systems is the most time-consuming part of
the SS method. We consider solving these linear systems by Block Krylov subspace methods tospeed up the SS method.
We solve the generalized eigenvalue problem derived from molecular orbital calculation of
dimerization of epidermal growth factor receptor by the SS method. The size n of the matrix is
96,234. The number L of right-hand sides of linear systems that appear in the SS method is 32.
In this calculation, 16 linear systems are solved. These linear systems are solved in parallelbecause these systems are mutually independent. As solvers for solving linear systems, the
COCG method and the Block COCG method with residual orthonormalization (Block
COCGrQ) are used. Numerical experiments are carried out on eight nodes of T2K-Tsukuba.
Each node consists of four quad-core AMD Opteron 8356 2.3GHz. Two sockets are assigned tosolve a linear system.
Fig. 8 shows the computation time of the SS method. The left graph shows the computation
time of the SS method when the COCG method is used. From Fig. 8 (a), we can see that thecomputation time for solving each linear system is very different. In this case, it takes about
7,300 seconds for computing eigenpairs. On the other hand, Fig. 8 (b) shows that the18

computation time for solving each linear system is almost the same when the Block COCGrQ
method is used. In this case, it takes about 900 seconds for computing eigenpairs. By using theBlock Krylov subspace methods, the computation time of the SS method can be reduced.

(a) COCG method.

(b) Block COCGrQ method

Fig. 8. Computation time of the SS method. ■: Matrix distribution, ■: Preconditioning, ■:Solving linear systems, ■: Extracting eigenpairs.

8. Conclusion
The following topics were investigated in the present study for the purposes of performancetuning.
- Parallel three-dimensional FFT algorithm for next-generation supercomputers- Performance evaluation of quadruple- and octuple-precision BLAS implementation on
GPUs- Acceleration of sparse matrix-vector multiplication by automatic storage format selection
with GPU- Performance evaluation of parallel two-dimensional FFT using Intel AVX instructions- Study of triple- and quadruple-precision floating point operations with GPU- Optimization and performance tuning technique for a numerical software on Multi-{core,socket, GPU} computing environments- Development of efficient Block Krylov subspace methods for linear systems with multiple
right-hand sides and their applications
Further research will be carried out into these topics in the future, in collaboration withcomputational material science researchers, in order to accelerate actual applications.

9. Publications
1)

2)

Y. Kubota and D. Takahashi, “Optimization of Sparse Matrix-Vector Multiplication by Auto
Selecting Storage Schemes on GPU”, Proc. 11th International Conference on Computational
Science and Its Applications (ICCSA 2011), Part II, Lecture Notes in Computer Science, No.
6783, pp. 547–561, Springer-Verlag, 2011.
T. Nakayama and D. Takahashi, “Implementation of Multiple-Precision Floating-Point

19

3)
4)

5)

Arithmetic Library for GPU Computing”, Proc. 23rd IASTED International Conference on
Parallel and Distributed Computing and Systems (PDCS 2011), pp. 343–349, 2011.
T. Otaki, K. Fujiyama, T. Imamura, S. Yamada and M. Machida, “Performance evaluation of
Eigen_sg+ASPEN.K2”, In the proceedings of HPCS2012, poster presentation, 2012.
D. Mukunoki and D. Takahashi, “Performance Comparison of Double, Triple and Quadruple
Precision Real and Complex BLAS Subroutines on GPUs”, Proc. ATIP/A*CRC Workshop
on Accelerator Technologies for High-Performance Computing: Does Asia Lead the Way?
(ATIP/A*CRC Workshop '12), pp. 788-790, 2012.
D. Takahashi, “An Implementation of Parallel 2-D FFT Using Intel AVX Instructions on
Multi-Core Processors”, Proc. 12th International Conference on Algorithms andArchitectures for Parallel Processing (ICA3PP 2012), Part II, Lecture Notes in Computer
Science, No. 7440, pp. 197-205, Springer-Verlag, 2012.

6)

Y. Nakamura, K.-I. Ishikawa, Y. Kuramashi, T. Sakurai, and H. Tadano, “Modified Block
BiCGSTAB for Lattice QCD”, Comput. Phys. Comm., Vol. 183, pp. 34—37, 2012.

7)

M. Naito, H. Tadano, and T. Sakurai, “A modified Block IDR(s) method for computing highaccuracy solutions”, JSIAM Letters, Vol. 4, pp. 25—28, 2012.
T. Imamura, S. Yamada, and M. Machida, “A High Performance SYMV Kernel on a
Fermi-core GPU”, High Performance Computing for Computational Science - VECPAR
2012, Lecture Note in Computer Science (LNCS) 7851, pp.59-71, 2013.
T. Imamura, T. Utsumi, X. Lin，S. Yamamda, and M. Machida, “Performance tuning for
SYMV kernels on multi-generation GPU’s, Fermi and Kepler”，IPSJ SIG-HPC，Vol.
2013-HPC-138(7), pp.1-7, 2013. (in Japanese)8)9)10) I. Yamazaki, H. Tadano, T. Sakurai, and T. Ikegami, “Performance comparison of parallel
eigensolver based on a contour integral method and a Lanczos method”, Parallel Computing,Vol. 39, Iss. 6—7, pp. 280—290, 2013.
11) T. Sakurai, Y. Futamura, and H. Tadano, “Efficient parameter estimation and implementationof a contour integral-based eigensolver”, J. Alg. Comput. Tech., Vol. 7, No. 3, pp. 249—269,
2013.
12) I. Yamazaki, A. Imakura, H. Tadano, and T. Sakurai, “A parameter estimation for thepreconditioning using space direct solvers for Krylov subspace methods based on theminimal residual approach”, Transactions of the Japan Society for Industrial and Applied
Mathematics, Vol. 23, No. 3, pp. 381—404, 2013 (in Japanese).
13) D. Takahashi, “Implementation of Parallel 1-D FFT on GPU Clusters”, Proc. 2013 IEEE 16th
International Conference on Computational Science and Engineering (CSE 2013), pp.
174—180, 2013.
14) T. Imamura, S. Yamada and M. Machida, "Eigen-G: GPU-based eigenvalue solver for
real-symmetric dense matrices", 10th International Conference on Parallel Processing andApplied Mathematics (PPAM2013), Lecture Note in Computer Science (LNCS) 8384, pp.
673--682, 2014.
15) A. Imakura, L. Du, and H. Tadano, “A weighted Block GMRES method for solving linear
systems with multiple right hand sides”, JSIAM Letters, Vol. 6, pp. 65-68, 2014.
16) S. Saito, H. Tadano, and A. Imakura, “Development of the Block BiCGSTAB(l) method forsolving linear systems with multiple right hand sides”, JSIAM Letters, Vol. 6, pp. 65—68,
2014.
17) H. Tadano, Y. Ishikawa, and A. Imakura, “Improvement of the accuracy of the approximate
solution of the Block BiCR method”, JSIAM Letters, Vol. 6, pp. 61—64, 2014.
18) T. Imamura, D. Mukunoki, S. Yamada, and M. Machida, “Implementation and evaluation of
CUDA-xSYMV”, IPSJ-SIHPC, Vol. 2014-HPC-146(14), pp.1-12, 2014. (in Japanese)19) T. Imamura, D. Mukunoki, S. Yamada and M. Machida, “Performance evaluation of the20

fastest GPU-based eigensolver by selecting CUDA-BLAS parameters”, IPSJ-SIGHPC, Vol.
2015-HPC-148(5), pp.1-9, 2015. (in Japanese)21

A01-3: Fast solution for ultra-large scale systems as a basisof computational materials science
S.-L. Zhang1, T. Sogabe2, I. Imahori1, T. Miyata1, Y. Yamamoto3
Nagoya University, 2Aichi Prefectural University,
3
The University of Electro-Communications

1

Finding novel complex correlation phenomena and clarifying the non-equilibrium
dynamics are of prime importance in the field of materials design. This researchgroup tackles the challenging problems from the viewpoints of numerical linear
algebra, optimization and high-performance computing.
The major purpose of the researches is to develop robust and efficient numericalalgorithms for solving large linear systems and eigenvalue problems in order toshed light on a breakthrough toward the challenging problems. For achieving the
purpose, the research group is organized by sophisticated experts as follows:
Professor Shao-Liang Zhang as the principal investigator; Professor Yusaku
Yamamoto as a project member in the field of high-performance computing forlarge eigenvalue problems; Associate professor Shinji Imahori as a project
member in the field of optimization; Associate professor Tomohiro Sogabe as aproject member in the field of numerical linear algebra for large linear systems;
Assistant professor Takafumi Miyata as a project member in the field of
numerical linear algebra for large eigenvalue problems.

1. Introduction
Computational science is now widely known to be a third mode of science for understanding
and predicting the behavior of artificial and natural systems with multiple degrees of freedom,
and it goes without saying that recent development of computational science highly depends oninnovation and breakthrough in the field of numerical linear algebra such as numericalalgorithms for linear systems and eigenvalue problems. In fact, millions of systems frequently
arise in the field of materials design through computics, and obtaining numerical solutions are
known to be one of the most time-consuming parts in the whole computation. Innovation andbreakthrough in terms of numerical linear algebra, therefore, will play an crucial role in further
development of materials design.
The primary purpose of the research is to develop high-performance numerical algorithms inorder to support disclosing hidden nature of the systems with multiple degrees of freedom. To
be more specific, the faster the better; the more accurate the better; and enhancing numericalmathematics for turning difficult problems into much simpler ones.

2. Novel high-performance numerical algorithms
In this section, the results on linear solvers, eigenvalue solvers and optimization solvers arereported.

22

2.1
2.1.1

Results on solvers for large linear systems
Subspace methods“Krylov subspace methods” is known as a Top10 algorithm of the 20th century [SIAM News,
Volume 33, Number 4], which is a iterative method for solving large and sparse linear systems.
The BiCG (Bi-Conjugate Gradient)-type method and the GMRES (Generalized Mimimal
REsidual) method are well known Krylov subspace methods for solving nonsymmetric linear
systems.
First, the result based on the GMRES method is described. As mentioned above, the GMRES
method is a representative Krylov solver. The method is proved to be optimal in terms of thenumber of iteration steps until convergence. It, however, is not practical in use due to growing
memory requirement and computational costs. Thus, restarted variant is widely used in practice.
Here, the restart means that after reaching at a prescribed iteration number m, the GMRES
method is restarted by setting the initial guess as the approximate solution at m iteration steps.
The restarted GMRES method, known as GMRES(m), is indeed useful in practice but there is a
disadvantage over the GMRES method in that the required number of iterations tends to glow
and at worst it may suffer from stagnation. In order to overcome the disadvantage, we focusedon the initial guess at every restart. The conventional GMRES(m) method uses the initial guess
to be the approximate solution obtained at the restart. We developed this part, i.e., we proposed
to set the initial guess to be the sum of the approximate solution and a modified vector. It is easy
to see that this approach includes the conventional one. Since there is a freedom in choosing themodified vector, we further investigated this effect in terms of matrix polynomials, and thenobtained how to choose the modified vector. As a result, we succeeded in obtaining themodified Look-Back GMRES(m) method. For some applications, the method outperformed the
GMRES(m) method.
Second, the result based on the BiCG method is described. The BiCG method is also a
representative Krylov solver and it is a cost-efficient algorithm in that the computational costs ateach iteration step is relatively small. The algorithm, unlike the GMRES(m) method, shows
irregular convergence behavior. To remedy the behavior, one of the most successful classes is
product-type Krylov subspace methods such as the BiCGStab (BiCG Stabilized) method andthe GPBiCG (Generalized Product-type BiCG) method. The key idea of the class is applying asuitable matrix polynomials to the BiCG residual vectors. The BiCGStab method uses a matrix
polynomial of degree 1 such that the resulting residual norm is locally minimized, and the
GPBiCG method uses a matrix polynomial of degree 2 for further stabilizing the BiCGStab
method. Recently, Sonneveld and van Gijzen proposed the Induced Dimension Reduction
method, so called IDR(s). The IDR(s) method finds approximate solution over Sonneveld
subspaces that include the subspace used by the BiCGStab method. In many applications the
IDR(s) method is known to be attractive. Recently, for further analyzing the property of theIDR(s) method, Sleijpen used a BiCG polynomial that is close to the IDR polynomial forclarifying the relation between the IDR(s) method and the BiCGStab method. Based on
Sleijpen’s analysis, we developed a product type Krylov subspace methods based on the BiCG
polynomial close to the IDR one, and then we obtained a modified BiCGStab method and amodified GPBiCG method. After the modification, in order to obtain father stabilization, we
proposed a Quasi-Minimal Residual variant of the IDRstab method. In fact, the variant often
gives smoother convergence behavior than the conventional one. In addition, we proposed amodified IDRstab method in order to reduce the number of DAXPY operations: addition of ascalar multiple of a real or complex vector to another real or complex vector using double

23

precision arithmetic.

2.1.2

Preconditioning techniquePreconditioning techniques are commonly used for boosting the speed of convergence of
Krylov subspace methods. The standard preconditioning techniques are transforming theoriginal system into another suitable one by means of multiplying the original systems by apreconditioning matrix. The resulting algorithm has a specific operation that corresponds tosolving the extra linear system with the preconditioning matrix at each iteration step. Since thepreconditioning matrix can be changed at each iteration step, there are some flexible
preconditioned Krylov solvers such as the FGMRES method, the GMRESR method, and the
GCR method with a variable preconditioning using the SOR (Successive Over-Relaxation)
method. When one uses the SOR method for the variable preconditioning, one cannot use theconventional SOR method for the case where the coefficient matrix has zero diagonal entries. In
order to remedy the difficulty, we modified a way of splitting the coefficient matrix andconstructed a variant of the generalized SOR method. From numerical experiments, we showed
the effectiveness of the variable preconditioning based on the variant of the generalized SOR
method. Furthermore, since solving the extra linear system is one of the most time-consuming
parts, we proposed a stopping criterion for solving the extra linear system from the viewpointsof the accuracy of the dot products arising in computing the coefficients of residual
polynomials.

2.1.3

Subspace methods for solving generalized shifted linear systems

Shifted linear systems with complex symmetric matrices arise in condensed matter physics
and the coefficient matrices are available in program package of ELSES (Extra Large Scale
Electronic Structure Calculation) led by Prof. Fujiwara at the university of Tokyo. The shifted
COCG method [1] is known as an efficient solver for the shifted linear systems.
We developed subspace methods for solving a generalized shifted linear system with complex
symmetric, which also arises in condensed matter physic. The key idea of the method is asfollows: first, transforming the generalized shifted linear systems into the standard shifted linear
systems; second applying the shifted COCG method with a bilinear form to the transformed
linear systems. The key idea made it possible to use the theory constructed in [1] for solving thegeneralized shifted ones.
For a numerical example, we use a matrix Au864 from large scale electronic structure
calculation, where the size of the coefficient matrix is 7776. The generalized shifted linear
systems to be solved is as follows(σ k S − H)x (k) = e 1 , k = 1,2,…,mwhere，σ k = −0.4+(k−1+i)/1000, m=1001, and H comes from Au864 and S is a given overlap
matrix. The result of comparison among the COCG method and our algorithms is shown in
Table.1.
Table1: Comparison among the COCG method and proposed methods ([2]).
CPU Time [s]SolverSeed system

Shifted system

Total

GSCOCG

486.2050.10

536.30

GSQMR

324.2261.61

385.83

GSQMRB

298.1649.28

347.44------

16867.8

COCG

24

We see from Table 1 that our algorithms outperform the COCG method in terms of the CPU
time. In fact, our algorithms are at least 30 times faster than the COCG method.

References:[1] T. Sogabe, T. Hoshi, S.-L. Zhang and T. Fujiwara, On a weighted quasi-residual minimization
strategy for solving complex symmetric shifted linear systems, Elec. Trans. Numer. Anal., 31(2008), pp. 126-140.
[2] T. Sogabe, T. Hoshi, S.-L. Zhang, and T. Fujiwara, Solution of generalized shifted linear systems
with complex symmetric matrices , J. Comput. Phys., 231 (2012), pp. 5669-5684.

2.1.4 On an improvement of Ye’s hybrid method for solving sphere constrained
quadratic opitimization
The sphere constrained quadratic optimization problem:

has a basic and important problem in optimization. Ye’s hybrid method [1] is one of the mostwell-known solvers. Ye’s hybrid method requires to solve the following shifted linear systemsof the form:

where Q is a symmetric positive definite matrix and the parameters are real. The fast solution of
the shifted linear systems is a bottleneck in the whole computation.
We proposed the implicit evaluation method for solving the above linear systems by using
shift-invariance property of Krylov subspaces. Let Q be a matrix bcsstk09 obtained from Matrix
Market (url: http://math.nist.gov/MatrixMarket/). Then the numerical result is given in Fig.1.

Fig.1. Comparison of explicit and implicit evaluations.
Step i (i=1,2,…,9) in Fig.1 means the computational costs for solving ith linear system to besolved. We see from Fig. 1 that total CPU time for the implicit evaluation method is almost thesame as the CPU time of step 1, which means 8 kinds of linear systems are efficiently solved.

25

Reference:[1] Y. Ye, A new complexity result on minimization of a quadratic function with a sphere constraint,
in C. Floudas and P. Pardalos eds., Recent Advances in Global Optimization, Princeton
University Press, NJ, 1992.

2.2

Subspace methods and high performance computing for eigenvalue problems
with applications

Eigenvalue problems as well as linear systems arise in rich variety of applications anddeveloping the efficient algorithms plays a crucial role in not only reducing the execution timeof the whole simulation but also deeper understanding physical phenomena. This section
describes the results of efficient algorithms suitable for computing specific eigenvalues and highperformance algorithm for all eigenvalues.

2.2.1

A fast solver for eigenvalue problems arising in photonic crystals

Photonic crystals have recently been intensively studied and the integrated circuits using
photons, instead of electrons, have received much attention. In the field of designing thenext-generation integrated circuits, computing eigenvalues located on the unit circle in thecomplex plane is required [1]. This requirement is much different from others. For example,
electromagneto dynamics and large scale electronic structure calculations require several
eigenvalues on the neighboring area of the targeted point in the complex plane, and some
important applications require extremal eigenvalues. Thus, algorithms for this purpose havebeen developed [2]. When it comes to the requirement in photonic crystals, it is a challenging
problem because all eigenvalues on the unit circle are required.
We focused on Sakurai-Sugiura (SS) method [3] as one of the most promising methods formeeting the requirement. The SS method finds all the eigenvalues inside the given region.
Based on this property, we extend the SS method to computing all the eigenvalues on the unit
circle. The key idea is to extend the SS method to computing all the eigenvalues on multiply
connected region. Then we set the multiply connected region to be double circles (inner andouter circles). By setting the two circles so that the unit circle is located between the two circles,
we can obtain all the eigenvalues on the region of the disc with hole. See Fig. 2 for details.
For numerical experiments, we used the eigenvalue problems from photonic crystals given by
Prof. Lu at Hong Kong University of Science and Technology, and showed that all eigenvalues
on the unit circles were obtained.
References:[1] Y. Huang, Y. Y. Lu, and S. Li, Analyzing photonic crystal waveguides by Dirichlet-to-Neumannmaps, J. Opt. Soc. Am. B, 24 (2007), pp. 2860-2867.
[2] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe and H. van der Vorst, eds., Templates for the solution of
Algebraic Eigenvalue Problems: A Practical Guide, SIAM, Philadelphia, 2000.
[3] T. Sakurai and H. Sugiura, A projection method for generalized eigenvalue problems using
numerical integration, J. Comput. Appl. Math., 159 (2003), pp. 119-128.

26

Fig. 2. Spectrum in the complex place. (Red circle: required eigenvalues.)

2.2.2

The Arnoldi (M,W,G) method for large eigenvalue problems

Large scale electronic structure calculations require solving large and sparse generalized
eigenvalue problems. For symmetric eigenvalue problems, the Lanczos method is widely known
and for nonsymmetric ones, the Arnoldi method is useful, which is an extension of the Lancos
methodWe succeeded in extending the Arnoldi method to solving generalized eigenvalue problems.
The resulting method is referred to as the Arnoldi(M, W, G) method. The relation among the theArnoldi(M, W, G) method, the Arnoldi method and the Lanczos method is as follows: Let thegeneralized eigenvalue problems be Ax=λBx. Then the Arnoldi(B, B*B, B) method corresponds
to the Arnoldi method and the Arnoldi(B, B, I) method corresponds to the Lanczos method.
From numerical experiments, we found the Arnoldi(I, I, I) method to be a method of choice fora problems arising in large scale electronic structure calculations

2.2.3

Eigenvalue problems arising in large scale electronic structure calculations

Recently, there is a strong need to compute the kth smallest eigenvalue for symmetriceigenvalue problems that arise in large scale electronic structure calculations.
One of the standard methods may be the Lancos method. In fact, the Lanczos method cancompute the smallest eigenvalue and neighboring eigenvalues. It is, however, difficult to tell
how many eigenvalues there are smaller than the obtained eigenvalue. This motivated us todesign an algorithm for avoiding the difficulty.
The key idea of the work is to use Sylvester's law of inertia together with a bisection search.
For a mall test problem, we successfully obtained the the kth smallest eigenvalue [1]．
Reference:[1] D. Lee, T. Miyata. T. Sogabe, T. Hoshi and S.-L. Zhang, Eigenvalue computation for a specificneed arising from electronic structure calculation, The 8th East Asia SIAM Conference, Taipei,
Taiwan, June 25-27, 2012.

27

2.2.4

Study on massively parallel algorithms for small scale symmetric eigenvalue
problemEigenvalue problem of a real symmetric matrix, Ax = λx, is an important problem thatfrequently arises in quantum chemistry and electronic structure calculations. We consider the
problem of computing all the eigenvalues and eigenvectors of a relatively small matrix (n ~10,000) on a massively parallel computer with thousands of computing nodes. The standard
algorithm for the symmetric eigenvalue problem is based on tri-diagonalization of the input
matrix, followed by the solution of the tridiagonal eigenvalue problem and back-transformationof the eigenvectors. While this approach is the best one among known algorithms in terms ofcomputational work, it is not efficient in our situation. This is because the tri-diagonalization
step is fine-grained and incurs many inter-processor communications. As a result, for relatively
small problems, the cost of inter-processor communication or synchronization becomesdominant and the speedup saturates around only a few hundred processors.
As an alternative to the tri-diagonalization based methods, we focus on the block-Jacobi
algorithm, which is a generalization of the well-known Jacobi algorithm for the symmetriceigenvalue problem. It works on a matrix partitioned into small blocks and tries to diagonalize
the matrix by eliminating the off-diagonal blocks by orthogonal transformations. Although thisalgorithm requires much more work than the tri-diagonalziation based algorithms, it has simpler
computational patterns and larger parallel granularity. Thanks to these features, the algorithm iseasier to optimize on modern high performance processors. Moreover, it is expected to attainbetter scalability for medium-size problems because the overhead of inter-processor
communication and synchronization is smaller. An additional feature of the block-Jacobi
algorithm is that most of the computation is done using the level-3 BLAS, or matrix-matrix
multiplication, which is a cache-friendly operation.
We implemented the block Jacobi method on the K computer and did some optimizations
such as hiding the inter-processor communications and tuning the solver of a 2 by 2 small
eigenproblem, which is the most time consuming part of the algorithm. According to ournumerical experiments, the resulting program shows good scalability on the K computer andoutperforms ScaLAPACK when solving an eigenvalue problem of size 10,000 using up to
10,000 computing nodes.
Another advantage of the block-Jacobi algorithm is that the convergence becomes faster
when the input matrix is nearly diagonal. In a certain type of first-principles molecular
dynamics method, the Hamiltonian matrix needs to be diagonalized at each time step. In thatcase, if we transform the Hamiltonian matrix at a time step with the (unitary) eigenvector matrix
at the preceding step, the resulting matrix becomes nearly diagonal. Hence, the convergence ofthe block-Jacobi method can be accelerated. We are planning to investigate this effect in moredetail and exploit it to speed up first-principles molecular dynamics methods.

2.2.5

Maximum Eigenvalue Problem and Local Search for Escher-Like Tilings

Maurits Cornelis Escher (1898-1972) is a Dutch graphic artist who is famous for impossible
structures, transformation prints and regular division of the plane. Inspired by his art works,
Kaplan and Salesin [1] introduced an optimization problem to design Escher-like tilings.
Koizumi and Sugihara [2] gave a new formulation of the problem and proposed an effective
algorithm for the problem; their method computes the maximum eigenvalue and corresponding
eigenvector to find an Escher-like tiling. A shortcoming of their method is sensitivity; for an
input line figure, different sets of points can approximate it, but they may output entirely28different tilings. In [2], good sets of points are generated manually through trial and error.
In this study [3], we propose two local search techniques to design good tilings automatically.
In the first method, a good set of points are automatically generated for an input line figure. For
the purpose, the features of the input figure are extracted and some fixed points are put on theboundary to keep the features. We then search good positions of remaining points forEscher-like tilings by local search technique. In the second method, constraint conditions arereconsidered. Koizumi and Sugihara uniquely determined a matrix to represent constraint
conditions for a tiling type. (It is known that there are 93 different tiling types for isohedraltilings.) Our method generates multiple matrices (i.e., many possibilities of shapes for tilings)
for each tiling type. By using these ideas, we design efficient algorithm for high-quality
Escher-like tilings.

Figure 1: Tiling of “Golden Dolphin”
We evaluate the proposed algorithm via computational experiments. In Figure 1 we show aresult of our method: (a) represents initial positions of points on the boundary of “Golden
Dolphin”, (b) is an optimal shape for (a) computed by the method in [2]. When our method
reached a local optimal solution, shape (c) that tiles the plane was output. Figure 1(d) is thetiling generated by this shape. From these figures, our method seems to work well for complex
input figures.
References:[1] C. S. Kaplan and D. H. Salesin, Escherization. Proceedings of SIGGRAPH, pp. 499-510 (2000).
[2] H. Koizumi and K. Sugihara, Maximum Eigenvalue Problem for Escherization. Graphs andCombinatorics, 27 (2011), pp. 431-439.
[3] S. Imahori and S. Sakai, A Local-Search Based Algorithm for the Escherization Problem. The
IEEE International Conference on Industrial Engineering and Engineering Management, (2012),
pp. 151-155.

3. Major publications
1)

M. E. A. El-Mikkawy and T. Sogabe, Notes on particular symmetric polynomials withapplications, Appl. Math. Comput., 215 (2010), pp. 3311-3317.

2)

M. E. A. El-Mikkawy and T. Sogabe, A new family of k-Fibonacci numbers, Appl. Math.
Comput., 215 (2010), pp. 4456-4461.

3)

T. Miyata, Y. Yamamoto and S.-L. Zhang, Performance Modeling of Multishift QR Algorithms
for the Parallel Solution of Symmetric Tridiagonal Eigenvalue Problems, (Proceeding) Lecture
Notes in Computer Science, 6082 (2010), pp. 401-412.

29

4)

Y. Mizuno, K. Ohi, T. Sogabe, Y. Yamamoto and Y. Kaneda, Four-point correlation functionof a passive scalar field in rapidly fluctuating turbulence: Numerical analysis of an exact
closure equation, Phys. Rev. E 82, 036316 (2010), pp. 1-9.

5)

A. Fukuda, Y. Yamamoto, M. Iwasaki, E. Ishiwata and Y. Nakamura, A Bäcklund t
ransformation between two integrable discrete hungry systems, Physics Letters A, 3753 (2011),
pp. 303-308.

6)

Y. Hirota, Y. Yamamoto and S.-L. Zhang, Cache optimization of a non-orthogonal joint
diagonalization method, JSIAM Letters, 3 (2011), pp. 9-12.
7) T. Miyata, Y. Yamamoto, T. Uneyama, Y. Nakamura and S.-L. Zhang, Optimization of theMultishift QR Algorithm with Coprocessors for Non-Hermitian Eigenvalue Problems, East
Asian Journal on Applied Mathematics, 1 (2011), pp. 187-196.
8) T. Sogabe and S.-L. Zhang, An extension of the COCR method to solving shifted linear
systems with complex symmetric matrices, East Asia J. on Appl. Math., 1 (2011), pp. 97-107.
9) Y. Tanaka, S. Imahori, M. Yagiura, Lagrangian-Based Column Generation for the Node
Capacitated In-Tree Packing Problem, Journal of the Operations Research Society of Japan, 54(2011), pp. 219-236.
10) H. Teng, T. Fujiwara, T. Hoshi, T. Sogabe, S.-L. Zhang and S. Yamamoto, Efficient andaccurate linear algebraic methods for large-scale electronic structure calculations with
non-orthogonal atomic orbitals, Phys. Rev. B 83, 165103 (2011), pp. 1-12.
11) Y. Yamamoto and Y. Hirota, A parallel algorithm for incremental orthogonalization based on
the compact WY representation, JSIAM Letters, 3 (2011), pp. 89-92.
12) A. Fukuda, Y. Yamamoto, M. Iwasaki, E. Ishiwata and Y. Nakamura, Error analysis for matrix
eigenvalue algorithm based on the discrete hungry Toda equation, Numerical Algorithms, 61(2012), pp. 243-260.
13) Y. Hama, A. Fukuda, Y. Yamamoto, M. Iwasaki, E. Ishiwata and Y. Nakamura, On Some
Properties of Discrete Integrable System Derived from Hungry Lotka-Volterra System ofMultiplicative Type, Journal of Math-for-Industry, 4(A) (2012), pp. 5-15.
14) T. Hoshi, S. Yamamoto, T. Fujiwara, T. Sogabe and S.-L. Zhang, J. Phys., An order-n electronic
structure theory with generalized eigenvalue equations and its application to a ten-million-atomsystem, J. Phys.: Condens. Matter, 24 (2012), pp. 1-5.
15) S. Imahori, Y. Miyamoto, H. Hashimoto, Y. Kobayashi, M. Sasaki and M. Yagiura, The
Complexity of the Node Capacitated In-Tree Packing Problem, Networks, 59 (2012), pp. 13-21.
16) S. Imahori and S. Sakai, A Local-Search Based Algorithm for the Escherization Problem, The
IEEE International Conference on Industrial Engineering and Engineering Management, (2012),
pp. 151-155.
17) A. Imakura, T. Sogabe and S.-L. Zhang, An Efficient Variant of the GMRES(m) Method Basedon the Error Equations, East Asian Journal on Applied Mathematics, 2 (2012), pp. 19-32.
18) D. Mori, Y. Yamamoto and S.-L. Zhang, Backward error analysis of the AllReduce algorithm
for householder QR decomposition, Japan Journal of Industrial and Applied Mathematics, 29(2012), pp. 111-130.
19) T. Sogabe, T. Hoshi, S.-L. Zhang and T. Fujiwara, Solution of generalized shifted linear systems
with complex symmetric matrices, Journal of Computational Physics, 231 (2012), pp.
5669-5684.
20) Y. Tanaka, S. Imahori, M. Sasaki and M. Yagiura, An LP-Based Heuristic Algorithm for theNode Capacitated In-Tree Packing Problem, Computers & Operations Research, 39 (2012),

30

pp. 637-646.
21) A. Fukuda, Y. Yamamoto, M. Iwasaki, E. Ishiwata and Y. Nakamura, On a shifted LR
transformation derived from the discrete hungry Toda equation, Monatshefte für Mathematik,
170 (2013), pp. 11-26.
22) A. Fukuda, E. Ishiwata, Y. Yamamoto, M. Iwasaki and Y. Nakamura, Integrable Discrete
Hungry Systems and Their Related Matrix Eigenvalues, Annali di Matematica Pura ed
Applicata, 192 (2013), pp. 423-445.
23) A. Fukuda, E. Ishiwata, Y. Yamamoto, M. Iwasaki and Y. Nakamura, Integrable discrete
hungry systems and their related matrix eigenvalues, Ann. Mat. Pura Appl., 192 (2013), pp423-445.
24) Y. Hu, H. Hashimoto, S. Imahori and M. Yagiura, Efficient construction heuristic algorithms forthe rectilinear block packing problem, Proceedings of the International Symposium on
Scheduling 2013, (2013), pp. 80-85.
25) A Imakura, T Sogabe and S.-L. Zhang, An Efficient Variant of the Restarted Shifted GMRES
Method for Solving Shifted Linear Systems, Journal of Mathematical Research with
Applications, 33 (2013), pp. 127-141.
26) J. Jia, T. Sogabe, and M.E.A. El-Mikkawy, Inversion of k-tridiagonal matrices with Toeplitz
structure, Comput. Math. Appl., 65 (2013), pp. 116-125.
27) D. Lee, T. Miyata, T. Sogabe, T. Hoshi and S.-L. Zhang, An interior eigenvalue problem fromelectronic structure calculations, Japan Journal of Industrial and Applied Mathematics, 30(2013), pp. 625- 633.
28) T. Miyata, S. Honda, R. Naito and S.-L. Zhang, A numerical approach to surface Green's
functions via generalized eigenvalue problems, Japan Journal of Industrial and Applied
Mathematics, 30 (2013), pp. 653- 660.
29) S. Imahori, Y. Chien, Y. Tanaka and M. Yagiura, Enumerating bottom-left stable positions forrectangle placements with qverlap, Journal of the Operations Research Society of Japan 57(2014), pp. 45-61.
30) Y. Yamamoto, Z. Lan and S. Kudo: Convergence analysis of the parallel classical block Jacobi
method for the symmetric eigenvalue problem, JSIAM Letters, 6 (2014), pp. 57-60.

5. Major invited talks
1)

S.-L. Zhang, A. Imakura and T. Sogabe, GMRES(m) Method with Look-Back-type Restart for
Solving Nonsymmetric Linear Systems, 7th International Conference on Scientific Computing
and Applications, Dalian, China, Jun. 13-15, 2010.

2)

S.-L. Zhang, T. Miyata, L. Du, T. Sogabe and Y. Yamamoto, A Projection Approach Based on
the Residue Theorem for Eigenvalues within a Multiply Connected Domain, The 2010 NIMS
Conference & The Third China-Japan-Korea Joint Conference on Numerical Mathematics,
Gangneung, Korea, Aug. 19-23, 2010.

3)

S.-L. Zhang and T. Miyata, A Numerical Method for Solving the Eigenvalues in a Multiple
Connected Domain, 2011 LHD Workshop, Beijing, China, Jan. 15, 2011.

4)

S.-L. Zhang and T. Miyata, A Projection Approach Based on the Residue Theorem for
Generalized Eigenproblems within a Multiply Connected Domain, Center of Mathematical31

Modeling and Scientific Computing, Chiao Tung University, Taiwan, Sep. 15, 2011.
5)

S.-L. Zhang and T. Sogabe, Iterative Solutions of Large Nonsymmetric Linear Systems ---Two
Attempts to Improve Bi-CG ---, Taida Institute for Mathematical Sciences National, National
Taiwan University, Taipei, Sep. 16, 2011.

6)

S.-L. Zhang, A. Imakura and T. Sogabe, Look-Back GMRES(m) for Solving NonsymmetricLinear Systems, Department of Mathematics, National Taiwan University, Taipei, Sep. 19, 2011.

7)

S.-L. Zhang,T. Miyata,T. Sogabe and T. Yamashita, An Arnoldi(M,W,G) Approach for
Generalized Eigenvalue Problems, National Center of Theoretical Science & Department ofMathematics, National Cheng Kung University, Tainan, Sep 20, 2011.

8)

S.-L. Zhang, A. Imakura and T. Sogabe, Look-Back GMRES(m) for Solving NonsymmetricLinear Systems, Department of Mathematics, National Sun Yat-sen University, Gaoxuong, Sep.
21, 2011.

9)

S.-L. Zhang,T. Miyata,T. Sogabe and T. Yamashita, An Arnoldi-like Approach for Generalized
Eigenvalue Problem, The 8th International Conference on NTmerical Optimization andNumerical Linear Algebra, Xiamen, China, Nov. 7-11, 2011.

10) S.-L. Zhang,T. Miyata and T. Sogabe, A Numerical Method for Eigenvalue Computations in a
Multiply Connected Domain, School of Mathematical Sciences, South China Normal University,
GuangZhou, China, Dec. 8, 2011.
11) S.-L. Zhang, A. Imakura and T. Sogabe, GMRES(m) method with Look-Back Restart for
Solving Nonsymmetric Linear Systems, the 2011 Workshop on Scientific Computing, Macau,
December 10-12, 2011.
12) S.-L. Zhang, A. Imakura and T. Sogabe, Look-Back GMRES(m) for Solving Large
Nonsymmetric Linear System, Numerical Linear Algebra ---Algorithms, Applications, andTraining, an NOW-JSPS joint seminar, Delft, Apr. 10-13, 2012.
13) T. Sogabe, T. Hoshi, S.-L. Zhang and T. Fujiwara, Numerical Methods for Generalized Shifted
Linear Systems with Complex Symmetric Matrices, 4th China-Japan-Korea Conference on
Numerical Mathematics, Otsu City, Shiga Prefecture, Japan, Aug. 25-28, 2012.
14) S.-L. Zhang, D. Lee, T. Miyata and T. Sogabe, An Intermediate Eigenvalue Problem inElectronic Structure Calculation, The 4th Workshop on Numerical Algebra and High
Performance Computation, Macau, Dec. 8-10, 2012.
15) Y. Yamamoto: A Nonlinear Eigenvalue Problem Arising in Theoretical Fluid Dynamics and Its
Solution Using Signed Singular Values, 2013 NCTS Workshop on Numerical Linear Algebraand High Performance Computing (2013 NLA-HPC), National Tsinghua University, Hsinchu,
Taiwan, Dec. 9-12, 2013.

32

Development of High Performance Parallel Eigensolvers
for the Needs in Computational Quantum Science

Tetsuya Sakurai
Department of Computer Science, University of Tsukuba, Tsukuba 305-8573
CREST, Japan Science and Technology Agency (JST), Kawaguchi 332-0012

1. Introduction
The large-scale eigenvalue problem is one of the common problems in computational design
and analysis of materials. High speed and high accuracy algorithms and their high performance
software are the fundamental technologies necessary to solve this problem efficiently. We haveinvestigated the complex moment-based eigensolver named the Sakurai-Sugiura method (the SS
method), which can selectively compute the eigenvalues in a certain region and thecorresponding eigenvectors based on a contour integral. Since the SS method has hierarchicalparallelism, it suits well for recent hierarchical architectures.
In order to efficiently apply the SS method to some applications in computational quantum
science on recent supercomputers, we specifically proceeded with 1) improving the algorithmand developing a highly parallel software implementation, and 2) developing interfaces andcomputational kernels specialized for each application.
Our Fortran-based software "z-Pares" supports the Reverse Communication Interface(RCI). Since RCI is independent of the particular distributed data structures of applications, RCI
makes it possible to provide a seamless connection in several applications. Our C Languagebased software is listed in the Scalable Library for Eigenvalue Problem Computation (SLEPc)
ver.3.4 as "CISS". In this research project, we also developed implementations using GPU for
speed-up of the software.

2. Achievements
The SS method selectively computes the eigenvalues in a certain region and the corresponding
eigenvectors. It constructs an invariant subspace associated with the target eigenvectors basedon a contour integral. The Arnoldi/Lanczos method, which is one of the most widely usedeigensolvers, constructs the subspace by sequential computation; on the other hand, the SS
method can construct the subspace by independent computation based on the contour integral.
Therefore, the SS method is expected to achieve higher parallel efficiency than traditional
eigensolvers.
The most time-consuming part of the SS method is to solve linear systems with multiple
right-hand sides on quadrature nodes. The SS method has the following three layers ofhierarchical parallelism; see Fig. 1.
Layer 1: Contour paths can be run independently.
Layer 2: Each linear system can be solved independently.
Layer 3: The linear systems can be solved in parallel.
A software implementation of this algorithmic hierarchy enables users to exploit theperformance of recent hierarchically structured architectures.
Based on cross-disciplinary cooperation with researchers in computational quantum
science, we have specifically proceeded as follows:33

Fig. 1: Hierarchical algorithm for solving sparse eigenproblems

Fig. 2: Parallel performance of an application of z-Pares to a matrix derived fromorder-N DFT.
(i) Improvement of the SS method and development of high performance software;(ii) Application to DFT code;(iii) Implementations using GPU for speed-up of the software.

2.1 Improvement of the SS method and development of high performance software
From research on applications in this application area, we fixed a specification of highperformance software z-Pares based on the SS method, and then developed z-Pares. z-Pares is
implemented in Fortran, which is widely used in computational physics, and used Message
Passing Interface (MPI) for implementations of distributed parallel processing. z-Pares supports
the Reverse Communication interface (RCI), which makes it possible to provide a seamlessconnection to several applications. z-Pares is available from http://zpares.cs.tsukuba.ac.jp/.

2.2 Application to DFT codes
2.2.1 Application to CONQUEST
We applied z-Pares to the order-N DFT code named CONQUEST. The test system is a SiGe hut
cluster system of 200,000 atoms. Single zeta basis is used as the basis set. Matrix size is34

778,292 and 223 eigenpairs around HOMO-LUMO are computed. We have executed theperformance evaluation on the supercomputer COMA at the University of Tsukuba. The total
number of nodes of COMA is 392. From these, we used up to 256 nodes during the application.
MUMPS is used as a linear solver. The parallel performance is shown in Fig. 2. The solid line
indicates parallel performance of z-Pares and the dotted line indicates an ideal scaling. From 4
nodes to 64 nodes, we used parallelism of quadrature points. In addition, from 64 nodes to 256,
we used parallelism of linear solver. We observe that up to 64 nodes parallel performance scales
almost linearly, whereas it degrades from 64 to 256 nodes. The linear scaling is caused by the
independence of computations of quadrature points. Parallel performance from 64 nodes isalmost equivalent to that of MUMPS. We observe that the parallel performance of the linear
solver is the key factor of highly parallel computational environment.

2.2.2 Application to RSDFTFor efficient parameter setting of the SS method, information of eigenvalue distribution (like theelectronic density-of-states) is needed. We have computed eigenvalue density of eigenproblems
derived from RSDFT code with contour integral based approach. To reduce the computational
cost of the approach, we proposed a variant of a method that efficiently computes block bilinear
form [10]. We have tested our approach on eigenproblems of SiNW with 9,924 atoms and107,292 atoms.

2.3 Implementations using GPU for speed-up of the software
Many of recent supercomputers listed in TOP500 have accelerators like GPU. In order tobenefit from this situation, we developed software for using GPU that is mainly based on zPares. The developed software for GPU uses the linear computational library Matrix Algebra on
GPU and Multicore Architectures (MAGMA) for some computation kernels of the SS method.
Fig. 3 shows the computational results of the developed software for computing 10%
eigenpairs of a 10,000 dimensional real symmetric generalized eigenvalue problem with random
dense matrices. We have executed the performance evaluation on the supercomputer HA-PACS
at University of Tsukuba. We used up to 16 GPUs of HA-PACS.
The most-time consuming part of the SS method is to solve independent linear systems: the"solve" and "fact" routines. Since these computations have trivial parallelism, total execution
time shows linear scaling specifically for the case of a small number of GPU.

Fig. 3: Calculation time of GPU implementation.

35

3.

Collaborators

Hiroto Tadano (University of Tsukuba) and Akira Imakura (University of Tsukuba)

4.

List of achievements

Paper[1] Akira Imakura, Lei DU, Tetsuya Sakurai, Error bounds of Rayleigh-Ritz type contour
integral-based eigensolver for solving generalized eigenvalue problems, Numerical
Algorithms, (accepted)[2] Lijiong Su, Akira Imakura, Hiroto Tadano and Tetsuya Sakurai, Improving theconvergence behaviour of the BiCGSTAB method by applying D-norm minimization,
JSIAM Letters, (accepted)[3] Shusaku Saito, Hiroto Tadano, Akira Imakura, Development of the Block BiCGSTAB(ell)method for solving linear systems with multiple right hand sides, JSIAM Letters, Vol. 6,
2014, pp. 65-68[4] Akira Imakura, Lei Du, Tetsuya Sakurai, A block Arnoldi-type contour integral spectral
projection method for solving generalized eigenvalue problems, Applied Mathematics
Letters, Vol. 32, 2014, pp.22-27[5] Hiroto Tadano, Youichi Ishikawa, Akira Imakura, Improvement of the accuracy of theapproximate solution of the Block BiCR method, JSIAM Letters, Vol. 6, 2014, pp. 61-64[6] Ikuro Yamazaki, Akira Imakura, Hiroto Tadano, Tetsuya Sakurai, A parameter estimation
for the preconditioning using sparse direct solvers for Krylov subspace methods based on
the minimal residual approach, Transactions of JSIAM, Vol. 23, 2013, pp.381-404 (inJapanese)[7] Akira Imakura, Tomohiro Sogabe and Shao-Liang Zhang, An efficient variant of therestarted shifted GMRES method for solving shifted linear systems, Journal of
Mathematical Research with Applications, Vol. 33, 2013, pp.127-141[8] Yuki Nagai, Yasushi Shinohara, Yasunori Futamura, Yukihiro Ota and Tetsuya Sakurai,
Numerical construction of a low-energy effective Hamiltonian in a self-consistent
Bogoliubov-de Gennes approach of superconductivity, J. Phys. Soc. Jpn., Vol. 82, 2013,
094701[9] T. Yano, Y. Futamura, T. Sakurai, Multi-GPU scalable implementation of a contourintegral-based eigensolver for real symmetric dense generalized eigenvalue problems,Proceedings of 8th International Conference on P2P, Parallel, Grid, Cloud and Internet
Computing (3PGCIC-2013), 2013, pp.121-127[10] Lei Du, Yasunori Futamura, Tetsuya Sakurai, Block conjugate gradient type methods forthe approximation of bilinear form, Computers & Mathematics with Applications, Vol. 66(12), 2014, pp.2446-2455
Award[11] Best paper award of 8th International Conference on P2P, Parallel, Grid, Cloud andInternet Computing (3PGCIC-2013), T. Yano, Y. Futamura, T. Sakurai, 201336

Hybrid mathematical principles
for ultra-large-scale massively parallel electronic structure calculations
Takeo Hoshi,

Tottori University

Abstract
Hybrid mathematical principles were introduced to large-scale electronic structure theories andrealized one-hundred-million atom (100-nm-scale) calculations on the K computer. Novellinear-algebraic theories were constructed as Krylov-subspace solvers for generalized shifted
linear equations ((zS-H)x=b) and were implemented in our order-N calculation code ‘ELSES’(http://www.elses.jp/) with modeled (tight-binding-form) systems. A high parallel efficiency
was obtained with up to the full core calculations on the K computer. The linear-algebraic
theory is general and, recently, was applied to ab initio quantum transport calculations on thereal space grid formalism. Moreover, hybrid mathematical priciples are constructed on iInternal
eigen-pair calculation with Sylvester's theorem of inertia and hybrid direct solver for
generalized eigen-value problem. Several application studies were carried out: The nano domainof sp2-sp3 nano-composite carbon solid was calculated for the nano-polycrystalline diamond, aultrahard material and typical domain shapes and boundary structures were obtained. An ionicliquid (N-Methyl-N-propylpiperidinium bis trifluoromethanesulfonyl imide), a battery-relatedmaterial, was calculated and the diffusion mechanism was revealed. Several amourphous-likestructures and bundles of conjugated polymers were calculated and characteristic localized πstates were obtained. The present work was carried out, with several members of Computics
Project: T. Sogabe, S.-L. Zhang, and T. Miyata (Nagoya U), T. Ono (Tsukuba U), asinter-disciplinary collaborations among physics, mathematics and the high-performancecomputation field or ‘Application-Algorithm-Architecture co-design'.

Fig. 1 (a) Concept of Application-Algorithm-Architecture co-design. (b) The
schematic workflow of the order-N algorithm based on the shifted Krylov subspace
theory.
1. Introduction
The current or next generation (exa-scale) computational physics requires‘Application-Algorithm-Architecture co-design', or inter-disciplinary collaborations among
physics, mathematics and the high-performance computation field. Figure 1(a) illustrates
our basic strategy for large-scale calculations with the co-design. The essential concept is tocombine new ideas among the three fields for a breakthrough. A schematic workflow is37

shown in Fig. 1(b). The mathematical solvers are general and can be shared, as middlewaresor ‘simulation kernels’, among applications not only in electronic structure calculations but
also in other computational science fields.

	  	 Fig. 2 (a) Benchmark for the order-N scaling. [7] The calculated system is an amorphous-likeconjugated polymer (aCP), poly-(9,9 dioctyl fluorene). (b) Benchmark for the parallel efficiency
or the strong scaling, on the K computer with one hundred million atoms. [1,2] The calculated
systems are aCP and an sp2-sp3 nano-composite carbon solid (NCCS).

2. Shifted Krylov subspace theory
One of our fundamental mathematical theories is novel linear algebraic algorithms
with Krylov subspace. The research was carried out with T. Sogabe and S.-L. Zhang(Nagoya U) within Computics Project. As shown in Fig. 1(b), the generalized shifted linear
equations ((zS-H)x=b) are solved instead of conventional generalized eigen-value problems.
A Krylov subspace is defined as the linear (Hilbert) space of Kn(A;b) ≡ span{ b, Ab,
A2b, ...., An-1b } with a given integer n, a square matrix A and a vector b. Krylov subspace isa general mathematical concept and gives the foundation of various iterative
linear-algebraic solvers, such as the conjugate gradient method.
Two ‘families’ of the solvers were constructed for the generalized shifted linear
equations. One is the solvers based on the collinear residual theorem (*1) and the other isthose based on the Galerkin principle. The solvers of the former family are the generalized
shifted Conjugate-Orthogonal Conjugate-Gradient method and the generalized shifted
Quasi-Residual Minimization method [6, 9]. The solvers of the latter family are thegeneralized Lanczos method, [9] the generalized Arnoldi method [9], the Arnoldi (M,W,G)
method [8] and the multiple Arnoldi method [7]. These solvers are compared in Ref. [7]. In
conclusion, the efficiency of the solvers depends on the system and purpose and themultiple Arnoldi method was adopted in the molecular dynamics simulations [7].
Figure 2 shows the benchmarks of the multiple Arnoldi method. The order-N scaling isshown in Fig. 2(a), in which the computational time is order-N or proportional to thenumber of atoms N. [7] Figure 2(b) shows the parallel efficiency on the K computer with
one hundred million atoms. [1,2] Since a system with one hundred million atoms is a
(126nm)3 region of silicon crystal, we call the calculations 100-nm-scale calculations. The
MPI/OpenMP hybrid parallelism is used. The results are shown for the nano-composite
carbon solid (NCCS) system (See Sec. 5) with N = 103, 219, 200 [2] and anamourphous-like conjugated polmer system, poly-(9,9 dioctyl fluorine) (See Sec.7), with N= 102,238,848 [1].

38The high parallel efficiency stems not only from the fundamental mathematical theory
but also from several programing techniques that do not affect the numerical results but
affects the computational time. In particular, programing techniques for saving memory andcommunication costs [7] are crucial for the massively parallel computations. Moreover, theparallel file I/O technique is sometimes important [1,7]. The K computer contains a parallel
file I/O function, called 'rank directory function' and was used in the present benchmark
with split input files [1,7].
The fundamental mathematical theory is general and, recently, was applied to the abinitio quantum transport calculation with real space grid, (*2) because the shifted linear
equaions appears on the conductor part of the systems. The research was carried out in thecollaboration with T. Ono (Tsukuba U) and his coworkers within Computics Project. As amethodological point, the shifted Conjugate-Orthogonal Conjugate-Gradient method is usedin the three-term recurrence formalism (*3), instead of the original two-term recurrence
formalism (*4) (*5). The three-term recurrence formalism is mathematically equivalent to
the two-term recurrence formalism but can save the memory consumption, when one uses
the seed switching technique, a techniques for optimal calculation among various situations.
(*1) W. A. Frommer, ‘BiCGStab(l) for Families of Shifted Linear Systems.’ Computing 70,
87 (2003).
(*2) S. Iwase, T. Hoshi, T. Ono, submitted.
(*3) S. Yamamoto, T. Sogabe, T. Hoshi, S.-L. Zhang and T. Fujiwara, 'Shifted
conjugate-orthogonal conjugate-gradient method and its application to double orbital
extended hubbard model', J. Phys. Soc. Jpn. 77,114713, 8pp, (2008).
(*4) R. Takayama, T. Hoshi, T. Sogabe, S.-L. Zhang, and T. Fujiwara, 'Linear algebraic
calculation of the Green's function for large-scale electronic structure theory', Phys. Rev. B
73, 165108, pp.1-9 (2006)(*5) The two-term and three-term recurrence formalisms are compared in several
mathematical papers from the viewpoint of numerical rounding error. For example, M. H.
Gutknecht and Z. Strakos, Accuracy of Two Three-term and Three Two-term Recurrencesfor Krylov Space Solvers, SIAM. J. Matrix Anal. & Appl., 22(1), 213–229 (2000).

Fig. 3 Iterative procedure for obtaining an internal eigen-value [3]. The target eigen value is the
HOMO level of an amorphous-like conjugated polymer. The number of the eigen values in the
search range is plotted as the function of iteration number. When the number of the eigen values isreduced to one, the iterative procedure will stop.

3. Internal eigen-pair calculation with Sylvester’s theorem of inertia

39

The theory was also constructed for calculating internal eigen pairs, sets of eigen
values and eigen vectors, with Sylvester’s theorem of inertia [3]. The research was carriedout with S.-L. Zhang, and T. Miyata (Nagoya U) within Computics Project. A popular needamong electronic structure calculations is to calculate specific eigen pairs, such as HOMO,
LUMO levels or levels near the Fermi level. The theory enables us to calculating the target
eigen pair with Sylvester’s theorem of inertia. The range that contains the target eigen value
will be determined in an iterative procedure, as in Fig. 3. After that, the precise eigen value
and the eigen vector is determined in the Arnoldi(M,W,G) method. [8]. The theory wastested with amourphous-like conjugated polymer and a nano metal. [3]Fig. 4 Benchmark of the hybrid direct solver on the K computer. (a) Strong scaling benchmark foran electronic structure problem in the matrix size of M=430,080 with up to 10,000 nodes. (b) The
explanation of the hybrid algorithms A-H used in (a).

4. Hybrid direct solver for generalized eigen-value problem
Optimally hybrid direct solvers were constructed for massively parallel generalized
eigenvalue computations. [10](*6) Though the direct solver consumes a heavy O(M3)operation cost with the matrix size M, it gives all the eigen pair exactly and forms a strong
need among electronic structure calculations. The generalized eigenvalue computation
procedure is decomposed into the two subprocedures of (i) the reducer to the standard
eigenvalue problem (SEP) and (ii) the solver of SEP. A hybrid solver is constructed, when aroutine is chosen for each subprocedure from the three parallel solver libraries of
ScaLAPACK, ELPA (*7) and EigenExa (*8). The strong scaling benchmark was carriedout on the K computer and other supercomputers for electronic structure calculation
problems in the matrix sizes of M = 104 -106. The hybrid solvers with the two newer
libraries, ELPA and EigenExa, give better benchmark results than the conventional
ScaLAPACK library, as shown in Fig. 4. Detailed performance analysis implies that theCholesky decomposition in the reducer can be a bottleneck in the next (post-petascale)computer and gives a guideline for further improvement on the algorithm. The code willappear online as a general middleware in Fig. 1(a).
(*6) H. Imachi and T. Hoshi, submitted.
(*7) A. Marek, et al., ’The ELPA Library - Scalable parallel eigenvalue solutions forelectronic structure theory and computational science’, J. Phys. Condensed Matter 26,

40

213201 (2014); http://elpa.rzg.mpg.de/(*8) T. Imamura, 'The EigenExa Library - High performance & scalable direct eigensolver
for large-scale computational science,’ ISC 2014, Leipzig, Germany (2014).'

Fig. 5 Nano-domain analysis of sp2-sp3 nano-composite carbon solid (NCCS) visualized with the
COHP or πCOHP analysis. [4] The sp2 (graphite-like) and sp3 (diamond-like) domains arevisualized in (a), while only the sp2-like domains are visualized in (b). A closeup of a sp2-sp3
domain boundary is shown in (c). The atoms at a typical domain boundary are shown in (d).

5. Nano-composite carbon solids
The large-scale calculation method was applied to nano-domain analysis of sp2-sp3
nano-composite carbon solids [4] for the formation process of the nano-polycrystalline
diamond (NPD) (*9). The NPD is a novel ultrahard material and is obtained by directconversion sintering process from graphite under high pressure and high temperature andhas characteristic 10-nm-scale lamellar-like structures. The NPD is of industrial importance,
because of its extreme hardness and strength that exceed those of conventional single
crystals. Sumitomo Electric Industries. Ltd. began commercial production from 2012.
Our simulation is motivated by the investigation of possible precursor structures in theformation process of NPD and the structures should be nano-scale composites of sp2(graphite-like) and sp3 (diamond-like) domains. Figures 5 shows such a result of our
finite-temperature simulation with a periodic boundary condition. Here the sp2 domains andthe sp3 domains are distinguished by a novel analysis method called 	  ‘π-crystalline
orbital Hamiltonian population (COHP) method’. The method is an analysis method for
local bond energy based on the Green’s function and is a theoretical extension of theoriginal COHP method. (*10) The πCOHP analysis is used, because sp2 domains can be
distinguished from sp3 domains by the presence of π bonds. Figures 5(a) and (b) are
drawn by the bond visualization by the original COHP and the πCOHP, respectively. In
Fig. 5(a), a bond is drawn for an atom pair according to the COHP analysis, and the sp2 andsp3 domains are visualized. In Fig. 5(b), on the other hand, a bond is drawn for an atom pairaccording to the πCOHP analysis, and only the sp2 domains are visualized. The analysis
concludes that layered domains are sp2 domains and non-layered domains are sp3 domains.
Figure 5(c) is a close-up of a domain boundary between sp2 and sp3 domains, which showsthat a layered domain forms a graphite-like structure and a non layered domain forms a
diamond-like structure, as expected from the πCOHP analysis. A typical domain interface,41

schematically shown in Fig. 5(d), is analyzed for the nature of each bond, so as to clarify
how theπCOHP analysis works. The atom labeled ‘A’ in Fig.5(d) is focused on. The atom
has three bonds and the two of them are located within an sp2 domain and have πbonds,
while the rest one is ‘interface’ bond between sp2 and sp3 domains and does not have a πbond. Such a domain shape or domain boundary structure will play crucial role in thestrength of NPD.
(*9) T. Irifune, A. Kurio, A. Sakamoto, T. Inoue, and H. Sumiya, ‘Ultrahard polycrystalline
diamond from graphite’, Nature 421, 599 (2003).
(*10) R. Dronskowski and P. E. Bloechl, ‘Crystal orbital Hamilton populations (COHP).
energy-resolved visualization of chemical bonding in solids based on density-functional
calculations’, J. Phys. Chem. 97 (1993) 8617; http://www.cohp.de/	 Fig. 6 (a)-(c) Snapshots of the finite-temperature dynamics simulation of the ionic liquid
PP13–TFSI at t=(a) 0 ps, (b) 10 ps, (c) 20 ps. (d) The density of states at t=20ps.

6. Ionic liquid
An industry-university collaboration research was carried out for the ionic liquid
N-Methyl-N-propylpiperidinium bis trifluoromethanesulfonyl imide (PP13–TFSI), a
battery-related material [5] (*11). The research is one for long time-scale phenomena, rather
than for large system. A finite-temerature simulation of the ion pairs in the period of 100 ps
at 303 K is shown in Fig. 6 and reveals the characteristic diffusion mechanism, since a
spontaneous dissociation of the ion-pair into the cation and the anion is not observed andthe ions diffuse as the ion-pair.
(*11) A co-auther of Ref. [5] is H. Yamasaki (Toyota Motor Corp.).
7. Conjugated polymer
Figure 7 shows preliminary results of conjugate polymers. [1,7,10] The research ispartially carried out as an industry-university collaboration. (*12) Conjugated polymer is ofgreat industrial importance among (opt)electronics and medical applications. The structures
are not in ideal crystals and the researches for condensed molecules in 100 nano-meter scale
are crucial. The calculated polymers in Fig. 4 are poly-(9,9 dioctyl fluorine) (PFO) andpoly-(para-phenylene vinylene) (PPV). Localized π wavefunctions are obtained from thetwisting (dihedral) angles between the benzene-like rings, such as in Fig.4 (b), and will playa key role for device property.
Now we are developing a transport calculation code based on the wavepacket
dynamics method for conjugated polymers and other organic materials, since the method isapplicable both to ballistic and non-ballistic transport.
.

42

(*12) Several structure data of amourphous-like conjugated polymers were provided by M
Ishida (Sumitomo Chemical Co.,Ltd) [1,7]Fig. 7 Caculattion of conjugated polymers; (a) Polymer structures of poly-(9,9 dioctyl
fluorine) (PFO)(left) and poly-(para-phenylene vinylene) (PPV) (right). (b) Localized
π wavefunction that appears on amourphous-like PFO polymers. [2] The inset
shows a schematic figure of the wave function. (c) Calculated PPV bundle with thelength L≒40 nm. [10]8. Collaborators within the Computics area
Tomohiro Sogabe, Takafumi Miyata, Shao-Liang Zhang (Nagoya U), Tomoya Ono (TsukubaU).

9. List of papers and invited talks
Original papers[1] Takeo Hoshi, Tomohiro Sogabe, Takafumi Miyata, Dongjin Lee, Shao-Liang Zhang, Hiroto
Imachi, Yoshifumi Kawai, Yohei Akiyama, Keita Yamazaki, Seiya Yokoyama, ‘Novel linear
algebraic theory and one-hundred-million-atom quantum material simulations on the K
computer’, PoS 202, 065, 13.pp (2014).
[2] Takeo Hoshi, Keita Yamazaki, Yohei Akiyama, “Novel linear algebraic theory andone-hundred-million-atom electronic structure calculation on the K computer”, JPS Conf. Proc.
1, 016004, 4pp. (2014).
[3] Dongjin Lee, Takafumi Miyata, Tomohiro Sogabe, Takeo Hoshi, Shao-Liang Zhang, “Aninterior eigenvalue problem from electronic structure calculations”, Japan J. Indust. Appl. Math.
30, pp 625-633 (2013).
[4] Takeo Hoshi, Yohei Akiyama, Tatsunori Tanaka and Takahisa Ohno,	  “Ten-million-atomelectronic structure calculations on the K computer with a massively parallel order-N theory”,
J. Phys. Soc. Jpn. 82, 023710, 4pp (2013).
[5] Shinya Nishino, Takeo Fujiwara, Hisatsugu Yamasaki, Susumu Yamamoto, Takeo Hoshi,"Electronic structure calculations and quantum molecular dynamics simulations of the ionicliquid PP13-TFSI ", Solid State Ionics 225, 22–25 (2012)[6] Tomohiro Sogabe, Takeo Hoshi, Shao-Liang Zhang, Takeo Fujiwara, "Solution ofgeneralized shifted linear systems with complex symmetric matrices ", J. Comp. Phys. 231,
5669-5684 (2012)[7] Takeo Hoshi, Susumu Yamamoto, Takeo Fujiwara, Tomohiro Sogabe, Shao-Liang Zhang,"An order-N electronic structure theory with generalized eigenvalue equations and its

43

application to a ten-million-atom system ", J. Phys.: Condens. Matter 24, 165502, 5pp. (2012)[8] Tatsuya Yamashita, Takafumi Miyata, Tomohiro Sogabe, Takeo Fujiwara, Shao-Liang
Zhang, " An Arnoldi (M, W, G) Method for Generalized Eigenvalue Problems", Trans of
JSIAM 21, 241 (2011) [JAPANESE].
[9] Honghui Teng, Takeo Fujiwara, Takeo Hoshi, Tomohiro Sogabe, Shao-Liang Zhang,
Susumu Yamamoto, "Efficient and accurate linear algebraic methods for large-scale electronic
structure calculations with nonorthogonal atomic orbitals ", Phys. Rev. B 83, 165103, 12pp,(2011)

Invited talks on international conference or workshop[10] Takeo. Hoshi and Hiroto Imachi, ‘Krylov subspace theories and one-hundred-million-atomelectronic structure calculations on the K computer’, International Symposium on Computics
Quantum Simulation and Design (ISC-QSD-2014), Koshiba Hall, University of Tokyo1-3, Dec.
(2014).
[11] Takeo Hoshi, Tomohiro Sogabe, Takafumi Miyata, and Shao-Liang Zhang, “Novel linear
algebraic theory and one-hundred-million-atom quantum material simulations on the K
computer”, The First International Workshop on Computational Science and Engineering ,National Taiwan University Taipei, Taiwan, 14-17., Oct., (2013).
[12] Takeo Hoshi, "Krylov subspace theories and ultra-large-scale electronic state calculations
on the K computer ", The 9th East Asia SIAM Conference The 2nd Conference on Industrial
and Applied Mathematics (EASIAM-CIAM 2013) , Institut Teknologi Bandung, Bandung, West
Java, Indonesia, 18-20., June (2013).
[13] Takeo Hoshi, Susumu Yamamoto, Takeo Fujiwara, Tomohiro Sogabe, Shao-Liang Zhang,
Yohei. Akiyama, Takahisa. Ohno, 'Ten-million-atom electronic structure calculations with novel
linear-algebraic algorithm and the K computer, International Symposium on Computics:
Quantum Simulation and Design (ISC-QSD), Osaka University Hall, 11-13. Oct. (2012)[14] Takeo Hoshi, Tomohiro Sogabe, Susumu Yamamoto, Shinya Nishino, Takeo Fujiwara,
Shao-Liang Zhang, "Novel linear-algebraic theories in large-scale electronic structure
calculation and its application to nano-materials ", Dresden-Kobe Joint Workshop on Electronic
Simulations for Nanosystems, Kobe University, 7. Mar. (2012)
Invited talks on domestic conference or workshop：totally 9 talks ( 1 talk in 2011, 2 talks in 2012,
4 talks in 2013, 1 talk in 2014, 1 talk in 2015).

44