We present optimal parallel QR factorization algorithms with reduced communication overhead. QR factorization is widely applied to solve various problems in numerical linear algebra. Our focus is on problems involving dense tall-and-skinny matrices in large-scale parallel distributed memory systems. Reducing data communication is essential for achieving high performance in parallel algorithms because the communication cost is much greater than the computation cost. To date, several QR factorization algorithms have been optimized to reduce communication costs. This paper provides alternative parallel QR factorization algorithms based on the LU-CholeskyQR algorithm. Numerical experiment results demonstrated the accuracy and performance of the developed algorithms against benchmarks. The results indicate that the new algorithms are numerically stable even for ill-conditioned problems, and some of these algorithms are faster than other conventional algorithms.

Mixed precision Krylov solvers with the Jacobi preconditioner often show significant convergence degradation when the Jacobi preconditioner is computed in low precision such as FP16 and BF16. It is found that this convergence degradation is attributed to loss of diagonal dominance due to roundoff errors in data conversion. To resolve this issue, we propose a new data conversion method, which is designed to keep diagonal dominance of the original matrix data. The proposed method is tested by computing the Poisson equation using the conjugate gradient method, the general minimum residual method, and the biconjugate gradient stabilized method with the FP16/BF16 Jacobi preconditioner on NVIDIA V100 GPUs. Here, the new data conversion is implemented by switching the round-nearest, round-up, round-down, and round-towards-zero intrinsics in CUDA, and is called once before the main iteration. Therefore, the cost of the new data conversion is negligible. When the coefficients of matrix is continuously changed by scaling the linear system, the conventional data conversion based on the round-nearest intrinsic shows periodic changes of the convergence property depending on the difference of the roundoff errors between diagonal and off-diagonal coefficients. Here, the period and magnitude of the convergence degradation depend on the bit length of significand. On the other hand, the proposed data conversion method is shown to fully avoid the convergence degradation, and robust mixed precision computing is enabled for the Jacobi preconditioner without extra overheads.

Although IEEE 754-2008 binary128 (with a 15-bit exponent and 113-bit significand, i.e., quadruple-precision) is not currently implemented on x86 in hardware, software emulation is available on some compilers. However, the performance is significantly slower compared to the binary64 operation, which is supported natively in hardware. This study proposes a fast implementation of matrix multiplication on matrices stored in the binary128 format on x86 CPUs. The proposed implementation utilizes the Ozaki scheme, which is an accurate matrix multiplication algorithm proposed by Ozaki et al. in 2012. This scheme enables one to perform most computations using the binary64 matrix multiplication (the DGEMM routine in Basic Linear Algebra Subprograms (BLAS)); it can exploit the high-performance of highly-optimized vendor BLAS. Although the achievable performance depends on the input matrices (the inner-product dimension, the absolute range, and the significand bit length), the proposed implementation can achieve better performance and accuracy compared to naive matrix multiplication performed using the GCC’s binary128 emulation in many cases. In addition, we discuss GPU acceleration, performance on reduced precision inputs, an implementation based on binary32 matrix multiplication (SGEMM), application to memory-intensive operations, and the possibility of a distributed parallel implementation.

Our performance benchmark of HPL-AI on the supercomputer Fugaku was awarded in the 55th top500 at ISC20. The effective performance was 1.42 EFlop/s, and the world's first achievement to exceed the wall of exascale in a floating-point arithmetic benchmark. Due to the novelty of HPL-AI, there are few guidelines for large systems and several drawbacks to the large-scale benchmark. It is not enough to replace FP64 operations solely to those on FP32 or FP16. At the least, we need thoughtful numerical analysis for lower-precision arithmetic and introduction of optimization techniques on extensive computing such as on Fugaku. In the poster, we give some comments on the accuracy, implementation, performance improvement, and report on the Exa-scale benchmark on Fugaku.

This research investigates the development of an FPGA-based sound field rendering system with the FDTD method, in which wave equations are directly implemented by the reconfigurable hardware, spatial blocking and temporal parallelism are applied to alleviate external memory bandwidth bottleneck and speed up computation. Compared with the software simulation carried out on a desktop machine with a Xeon Gold 6212U processor and 512 GB DRAMs, the FPGA system achieves about 13 times speedup in computation performance.

The multi-scale full-f simulation of the next generation experimental fusion reactor ITER based on a five dimensional (5D) gyrokinetic model is one of the most computationally demanding problems in fusion science. In this work, a Gyrokinetic Toroidal 5D Eulerian code (GT5D) is accelerated by a new mixed-precision communication-avoiding (CA) Krylov method. The bottleneck of global collective communication on accelerated computing platforms is resolved using a CA Krylov method. In addition, a new FP16 preconditioner, which is designed using the new support for FP16 SIMD operations on A64FX, reduces both the number of iterations (halo data communication) and the computational cost. The performance of the proposed method for ITER size simulations with ~0.1 trillion grids on 1,440 CPUs/GPUs on Fugaku and Summit shows 2.8× and 1.9× speedups respectively from the conventional non-CA Krylov method, and excellent strong scaling is obtained up to 5,760 CPUs/GPUs.

Numerical weather prediction (NWP) supports our daily lives. Weather models require higher spatiotemporal resolutions to prepare for extreme weather disasters and reduce the uncertainty of predictions. The accuracy of the initial state of the weather simulation is also critical; thus, we need more advanced data assimilation (DA) technology. By combining resolution and ensemble size, we have achieved the world’s largest weather DA experiment using a global cloud-resolving model and an ensemble Kalman filter method. The number of grid points was $\sim$4.4 trillion, and 1.3 PiB of data was passed from the model simulation part to the DA part. We adopted a data-centric application design and approximate computing to speed up the overall system of DA. Our DA system, named NICAM-LETKF, scales to 131,072 nodes (6,291,456 cores) of the supercomputer Fugaku with a sustained performance of 29 PFLOPS and 79 PFLOPS for the simulation and DA parts, respectively.

We herein propose an efficient implementation of tridiagonalization (TRD) for small matrices on manycore CPUs. Tridiagonalization is a matrix decomposition that is used as a preprocessor for eigenvalue computations. Further, TRD for such small matrices appears even in the HPC environment as a subproblem of large computations.

Most computer programs are written using either or both 32- and 64-bit floating-point formats (FP32 and FP64). To save memory space, improve the computation speed, and improve energy efficiency by eliminating waste bit data in floating-point data, we need reduced-precision formats that can represent various and lower floating-point precisions. This paper proposes the implementation of reduced-precision floatingpoint formats, which provide different and shorter significand lengths compared with the IEEE formats on graphics processing units. This paper demonstrates the performance of some linear algebra kernels (AXPY, GEMV, and CSRMV) using reduced-precision formats on GeForce GTX 980. Our implementations with the reduced-precision formats outperformed the optimized vendor implementations for FP32 and FP64.

This paper proposes a fast implementation method for the general matrix-vector multiplication (GEMV) routine, which is one of the level-2 Basic Linear Algebra Subprograms (BLAS) subroutines, for a column-major and non-transposed matrix on NVIDIA Kepler architecture graphics processing units (GPUs). We began by implementing the GEMV kernel using typical blocking techniques for shared-memory and register along with 128-bit vector load/store instructions. In our initial investigation, we found that even though the kernel could approach actual peak GPU throughput at some matrix sizes, performance fluctuates periodically depending on the problem size. In our next step, we investigated the reason for the fluctuations using a performance model based on a thread-block scheduling mechanism, and then created a method of determining optimal thread-block sizes that avoids those fluctuations. As the results show, when run on two Kepler architecture GPUs, our single-precision GEMV (SGEMV) routine achieved better performance in terms of both throughput and performance stability (with respect to the problem size) when compared to existing implementations: CUBLAS 6.5, MAGMA 1.4.1 and KBLAS 1.0. Our implementation techniques can be used not only for SGEMV but also double-precision (DGEMV), single-complex (CGEMV), and double-complex (ZGEMV). While this paper discusses primarily Kepler architecture, we also explore the performance of proposal implementation on Maxwell architecture, which is the next generation of Kepler architecture.

One of the most fascinating issues in modern condensed matter physics is to understand highly-correlated electronic structures and propose their novel device designs toward the reduced carbon-dioxide future. Among various developed numerical approaches for highly-correlated electrons, the density matrix renormalization group (DMRG) has been widely accepted as the most promising numerical scheme compared to Monte Carlo and exact diagonalization in terms of accuracy and accessible system size. In fact, DMRG almost perfectly resolves one-dimensional chain like long quantum systems. In this paper, we suggest its extended approach toward higher-dimensional systems by high-performanc computing techniques. The computing target in DMRG is a huge non-uniform sparse matrix diagonalization. In order to efficiently parallelize the part, we implement communication step doubling together with reuse of the mid-point data between the doubled two steps to avoid severe bottleneck of all-to-all communications essential for the diagonalization. The technique is successful even for clusters composed of more than 1000 cores and offers a trustworthy exploration way for two-dimensional highly-correlated systems.

In order to study quantum many-body problems, we develop two matrix diagonalization codes, which solve only the ground state and all quantum states, respectively. The target model in both codes is the Hubbard model with confinement potential which describes an atomic Fermi gas loaded on an optical lattice and partly High-Tc cuprate superconductor. For the former code, we obtain 18.692TFlops (57% of the peak) as the best performance on the Earth Simulator when calculating the ground state of 100-billion dimensional matrix. From these large-scale calculations, we find atomic-scale inhomogeneous superfluid state which is now a challenging subject for physicists. For the latter code, we succeed in solving the matrix whose dimension is 375,000 with locally 24.6TFlops (75% of the peak). The calculations reveal that a change from Schrodinger's cat to classical like one can be controlled by tuning the interaction. This is a marked contrast to the general concept.

In order to study a possibility of superfluidity in trapped atomic Fermi gases loaded on optical lattices, we implement an exact diagonalization code for the trapped Hubbard model on the Earth Simulator. Comparing two diagonalization algorithms, we find that the performance of the preconditioned conjugate gradient (PCG) method is 1.5 times superior to the conventional Lanczos one since the PCG method can conceal the communication overhead much more efficiently. Consequently, the PCG method shows 16.447 TFlops (50.2% of the peak) on 512 nodes. On the other hand, we succeed in solving a 159-billion-dimensional matrix by using the conventional Lanczos method. To our knowledge, this dimension is a world-record. Numerical results reveal that an unconventional type of superfluidity specific to the confined system develops under repulsive interaction.