BLAS Kernels Research Articles

Hierarchical approach for deriving a reproducible unblocked LU factorization

We propose a reproducible variant of the unblocked LU factorization for graphics processor units (GPUs). For this purpose, we build upon Level-1/2 BLAS kernels that deliver correctly-rounded and reproducible results for the dot (inner) product, vector scaling, and the matrix-vector product. In addition, we draw a strategy to enhance the accuracy of the triangular solve via iterative refinement. Following a bottom-up approach, we finally construct a reproducible unblocked implementation of the LU factorization for GPUs, which accommodates partial pivoting for stability and can be eventually integrated in a high performance and stable algorithm for the (blocked) LU factorization.

Optimized sparse Cholesky factorization on hybrid multicore architectures

We present techniques for supernodal sparse Cholesky factorization on a hybrid multicore platform consisting of a multicore CPU and GPU. The techniques are the subtree algorithm, pipelining and multithreading. The subtree algorithm [15] minimizes PCIe transmissions by storing an entire branch of the elimination tree in the GPU memory (the elimination tree is a tree data structure describing the workflow of the factorization), and also reduces the total kernel launch time by launching BLAS kernels in batches. The pipelining technique overlaps the execution of GPU kernels and PCIe data transfers. The multithreading technique [17] creates multiple threads for both the CPU and the GPU, to utilize concurrency of the elimination tree. Our experimental results on a platform consisting of an Intel multicore processor along with an Nvidia GPU indicate a significant improvement in performance and energy over CHOLMOD (SuiteSparse 4.5.3), a sparse algorithm, after these techniques are applied.

Reordering Strategy for Blocking Optimization in Sparse Linear Solvers

Solving sparse linear systems is a problem that arises in many scientific applications, and sparse direct solvers are a time-consuming and key kernel for those applications and for more advanced solvers such as hybrid direct-iterative solvers. For this reason, optimizing their performance on modern architectures is critical. The preprocessing steps of sparse direct solvers---ordering and block-symbolic factorization---are two major steps that lead to a reduced amount of computation and memory and to a better task granularity to reach a good level of performance when using BLAS kernels. With the advent of GPUs, the granularity of the block computation has become more important than ever. In this paper, we present a reordering strategy that increases this block granularity. This strategy relies on block-symbolic factorization to refine the ordering produced by tools such as Metis or Scotch, but it does not impact the number of operations required to solve the problem. We integrate this algorithm in the PaSt...

Acceleration of GPU-based Krylov solvers via data transfer reduction

Krylov subspace iterative solvers are often the method of choice when solving large sparse linear systems. At the same time, hardware accelerators such as graphics processing units continue to offer significant floating point performance gains for matrix and vector computations through easy-to-use libraries of computational kernels. However, as these libraries are usually composed of a well optimized but limited set of linear algebra operations, applications that use them often fail to reduce certain data communications, and hence fail to leverage the full potential of the accelerator. In this paper, we target the acceleration of Krylov subspace iterative methods for graphics processing units, and in particular the Biconjugate Gradient Stabilized solver that significant improvement can be achieved by reformulating the method to reduce data-communications through application-specific kernels instead of using the generic BLAS kernels, e.g. as provided by NVIDIA’s cuBLAS library, and by designing a graphics processing unit specific sparse matrix-vector product kernel that is able to more efficiently use the graphics processing unit’s computing power. Furthermore, we derive a model estimating the performance improvement, and use experimental data to validate the expected runtime savings. Considering that the derived implementation achieves significantly higher performance, we assert that similar optimizations addressing algorithm structure, as well as sparse matrix-vector, are crucial for the subsequent development of high-performance graphics processing units accelerated Krylov subspace iterative methods.

A BLAS-3 Version of the QR Factorization with Column Pivoting

The QR factorization with column pivoting (QRP), originally suggested by Golub [Numer. Math., 7 (1965), 206--216], is a popular approach to computing rank-revealing factorizations. Using Level 1 BLAS, it was implemented in LINPACK, and, using Level 2 BLAS, in LAPACK. While the Level 2 BLAS version delivers superior performance in general, it may result in worse performance for large matrix sizes due to cache effects. We introduce a modification of the QRP algorithm which allows the use of Level 3 BLAS kernels while maintaining the numerical behavior of the LINPACK and LAPACK implementations. Experimental comparisons of this approach with the LINPACK and LAPACK implementations on IBM RS/6000, SGI R8000, and DEC AXP platforms show considerable performance improvements.

A block variant of the GMRES method on massively parallel processors

This paper presents a block variant of the GMRES method for solving general unsymmetric linear systems. This algorithm generates a transformed Hessenberg matrix by solely using block matrix operations and block data communications. It is shown that this algorithm with block size s, denoted by BVGMRES( s, m), is theoretically equivalent to the GMRES( s, m) method. The numerical results demonstrate that this algorithm can be more efficient than the standard GMRES method on a cache based single CPU computer with optimized BLAS kernels. Furthermore, the gain in efficiency is more significant on MPPs due to both efficient block operations and efficient block data communications. Preliminary numerical results on some real-world problems also show that this algorithm may be stable up to some reasonable block size.

A block varaint of the GMRES method for unsymmetric linear systems

Iterative methods that take advantage of efficient block operations and block communications are popular research topics in parallel computation. These methods are especially important on Massively Parallel Processors (MPP). This paper presents a block variant of the GMRES method for solving general unsymmetric linear systems. It is shown that the new algorithm with block sizes, denoted by BVGMRES (s.m), is theoretically equivalent to the GMRES (s·m) method. The numerical results show that this algorithm can be more efficient than the standard GMRES method on a cache besed single CPU computer with optimized BLAS kernels. Furthermore, the gain in efficiency is more significant on MPPs due to both efficient block operations and efficient block data communications. Our numerical results also show that in comparison to the standard GMRES method, the more PEs that are used on an MPP, the more efficient the BVGMRES(s,m) algorithm is.

Performance of level 3 BLAS kernels in a dynamically partitioned data-flow environment

The Dynamically Partitioned Data-Flow (DPDF) model is based on an original analysis concept of the data dependency graph at the instruction level. Instead of a breadth first analysis, as in a classical Data-Flow Model, we execute instructions along data-dependent paths. As a consequence, data locality can be exploited by reusing results between the execution of consecutive instructions. In addition, the different paths are not statically defined but arise from a dynamical partitioning of the graph. This model presents the advantage to support very small cost dynamic scheduling and multitasking strategies. In order to study the efficiency of this new model, a first architecture has been defined. This architecture is currently limited to a single processor with one serial processing unit but four graph analyzing units (called prefetch units). Each of these prefetch units is able to build dynamically its own execution path inside the Data-Flow graph of an application. The efficiency of this architecture is studied on a numerical benchmark composed of a subset of the Livermore loops and of three routines of the Level 3 BLAS (GEMM, SYRK and TRSM). Our goal in these experimentations is to demonstrate the ability of the four prefetch units to feed the ALU.

Highly nonnormal eigenproblems in the aeronautical industry

We consider a large-scale nonnormal eigenvalue problem that occurs in flutter analysis. Matrices arising in such problems are usually sparse, of large order, and highly nonnormal. We use the incomplete Arnoldi method associated with the Tchebycheff acceleration in order to compute a subset of the eigenvalues and their associated eigenvectors. This method has been parallelized using BLAS kernels and has been tested on various vector and parallel machines. We also studied the stability of the eigenproblem by perturbing the original matrix and verified that the spectral instability increases with the size of the problem. This work has been conducted at CERFACS in cooperation with the Aerospatiale Avions (Structural research and development department).

Comparisons of Gaussian elimination algorithms on a cray Y-MP

The results of various implementations of Gaussian elimination on full matrices on a single processor Cray Y-MP are presented and discussed. It is shown that when the manufacturer supplied BLAS kernels are used, the difference between thebest versions of level 2 BLAS and level 3 BLAS (blocked) implementations is almost negligible. It is suggested that to improve the performance of blocked Gaussian elimination it is possible to utilize Strassen's matrix multiplication algorithm.

Block-Cholesky for parallel processing

We concentrate on the Cholesky factorization of A = LL T , where A is a positive definite symmetric matrix and L is a lower triangular matrix. A blocked algorithm based on Level 3 BLAS is discussed. When using Level 3 BLAS kernels in a multiprocessing mode, one can parallelize within each kernel, or can obtain parallelism by performing different matrix-matrix operations on different processors. We apply parallelism over the blocks. We study the amount of parallelism and we discuss the data dependency graph. The SCHEDULE package is used to obtain a portable scheduling of the tasks. Numerical results of our method are presented and compared with the results for a block algorithm parallelized within the Level 3 BLAS kernels.

Use of Level 3 Blas in Lu Factorization in a Multiprocessing Environment On Three Vector Multiprocessors: the Alliant Fx/80, the Cray-2, and the Ibm 3090 Vf

We study various implementations of block Gaussian elimination on full matrices and examine their perfor mance on three parallel computers, the Alliant FX/80, the CRAY-2, and the IBM 3090-400/VF. These imple mentations are expressed in terms of Level 3 BLAS matrix-matrix kernels. We consider the use of parallel Level 3 BLAS kernels and compare the parallelism ob tained within the computational kernels with that ob tained when parallelizing over the kernels. We show that the use of parallel Level 3 BLAS allows portability without sacrifice of efficiency, even in a parallel envi ronment, and that high speeds can be obtained if tuned versions of the kernels are available.

Level 3 Blas in Lu Factorization On the Cray-2, Eta-10P, and Ibm 3090-200/Vf

We study various implementations of block Gaussian elimination on full matrices and examine their perfor mance on three vector supercomputers, the CRAY-2, the ETA-10P, and the IBM 3090-200/VF. We show that the use of Level 3 BLAS kernels allows portability without sacrifice of efficiency and that good speeds can be ob tained if tuned versions of the kernels are available. In deed our results show that without using any assembler language outside the kernels we can approach the per formance of assembler-coded routines on all machines.