BLAS Research Articles

BLIS: A Framework for Rapidly Instantiating BLAS Functionality

The BLAS-like Library Instantiation Software (BLIS) framework is a new infrastructure for rapidly instantiating Basic Linear Algebra Subprograms (BLAS) functionality. Its fundamental innovation is that virtually all computation within level-2 (matrix-vector) and level-3 (matrix-matrix) BLAS operations can be expressed and optimized in terms of very simple kernels. While others have had similar insights, BLIS reduces the necessary kernels to what we believe is the simplest set that still supports the high performance that the computational science community demands. Higher-level framework code is generalized and implemented in ISO C99 so that it can be reused and/or reparameterized for different operations (and different architectures) with little to no modification. Inserting high-performance kernels into the framework facilitates the immediate optimization of any BLAS-like operations which are cast in terms of these kernels, and thus the framework acts as a productivity multiplier. Users of BLAS-dependent applications are given a choice of using the traditional Fortran-77 BLAS interface, a generalized C interface, or any other higher level interface that builds upon this latter API. Preliminary performance of level-2 and level-3 operations is observed to be competitive with two mature open source libraries (OpenBLAS and ATLAS) as well as an established commercial product (Intel MKL).

A parallel generalized relaxation method for high-performance image segmentation on GPUs

Fast and scalable software modules for image segmentation are needed for modern high-throughput screening platforms in Computational Biology. Indeed, accurate segmentation is one of the main steps to be applied in a basic software pipeline aimed to extract accurate measurements from a large amount of images. Image segmentation is often formulated through a variational principle, where the solution is the minimum of a suitable functional, as in the case of the Ambrosio–Tortorelli model. Euler–Lagrange equations associated with the above model are a system of two coupled elliptic partial differential equations whose finite-difference discretization can be efficiently solved by a generalized relaxation method, such as Jacobi or Gauss–Seidel, corresponding to a first-order alternating minimization scheme. In this work we present a parallel software module for image segmentation based on the Parallel Sparse Basic Linear Algebra Subprograms (PSBLAS), a general-purpose library for parallel sparse matrix computations, using its Graphics Processing Unit (GPU) extensions that allow us to exploit in a simple and transparent way the performance capabilities of both multi-core CPUs and of many-core GPUs. We discuss performance results in terms of execution times and speed-up of the segmentation module running on GPU as well as on multi-core CPUs, in the analysis of 2D gray-scale images of mouse embryonic stem cells colonies coming from biological experiments.

CheMPS2: Improved DMRG-SCF routine and correlation functions

CHEMPS2, our spin-adapted implementation of the density matrix renormalization group (DMRG) for ab initio quantum chemistry (Wouters et al., 2014), has several new features. A speed-up of the augmented Hessian Newton-Raphson DMRG self-consistent field (DMRG-SCF) routine is achieved with the direct inversion of the iterative subspace (DIIS). For extended molecules, the active space orbitals can be localized by maximizing the Edmiston-Ruedenberg cost function. These localized orbitals can be ordered according to the topology of the molecule by approximately minimizing the bandwidth of the exchange matrix with the Fiedler vector. The electronic structure can be analyzed by means of the two-orbital mutual information, spin, spin-flip, density, and singlet diradical correlation functions.

Bayesian hierarchical models for spatially misaligned data in R

Summary Spatial misalignment occurs when at least one of multiple outcome variables is missing at an observed location. For spatial data, prediction of these missing observations should be informed by within location association among outcomes and by proximate locations where measurements were recorded. This study details and illustrates a Bayesian regression framework for modelling spatially misaligned multivariate data. Particular attention is paid to developing valid probability models capable of estimating parameter posterior distributions and propagating uncertainty through to outcomes' predictive distributions at locations where some or all of the outcomes are not observed. Models and associated software are presented for both Gaussian and non‐Gaussian outcomes. Model parameter and predictive inference within the proposed framework is illustrated using a synthetic and forest inventory data set. The proposed Markov chain Monte carlo samplers were written in c++ and leverage R's Foreign Language Interface to call fortran blas (Basic Linear Algebra Subprograms) and lapack (Linear Algebra Package) libraries for efficient matrix computations. The models are implemented in the spMisalignLM and spMisalignGLM functions within the spBayes r package available via the Comprehensive R Archive Network (cran) (http://cran.r-project.org).

Performance evaluation of kernel fusion BLAS routines on the GPU: iterative solvers as case study

Programmers usually implement iterative methods that solve partial differential equations by expressing them using a sequence of basic kernels from libraries optimized for the graphics processing unit (GPU). The global runtime of the resulting combination is often penalized by the smallest and most inefficient vector operations. To improve the GPU exploitation, we identify and analyze the potential kernels to be fused according to the data dependence, data type and size, and GPU resources. This paper provides an extensive analysis of the impact of fusing vector operations [level 1 of Basic Linear Algebra Subprograms (BLAS)] on the performance of the GPU. The experimental evaluation shows that this optimization provides noticeable improvement especially for kernels with lower memory requirements and on more modern GPUs. It is worth noting that the fused BLAS operations can be very useful to help programmers efficiently code iterative methods to solve large linear systems of equations for the GPU. Iterative methods such as biconjugate gradient method (BCG) are one of the examples that can benefit from this optimization strategy. Indeed, kernel fusion of vector routines makes the most efficient GPU implementation of BCG run between $$1.09\times $$ 1.09 × and $$1.27\times $$ 1.27 × faster on three GPUs of different characteristics.

Code Generation and Optimization of Distributed-memory Dense Linear Algebra Kernels

Design by Transformation (DxT) is an approach to software development that encodes domain-specific programs as graphs and expert design knowledge as graph transformations. The goal of DxT is to mechanize the generation of highly-optimized code. This paper demonstrates how DxT can be used to transform sequential specifications of an important set of Dense Linear Algebra (DLA) kernels, the level-3 Basic Linear Algebra Subprograms (BLAS3), into high-performing library routines targeting distributed-memory (cluster) architectures. Getting good BLAS3 performance for such platforms requires deep domain knowledge, so their implementations are manually coded by experts. Unfortunately, there are few such experts and developing the full variety of BLAS3 implementations takes a lot of repetitive effort. A prototype tool, DxTer, automates this tedious task. We explain how we build on previous work to represent loops and multiple loop-based algorithms in DxTer. Performance results on a BlueGene/P parallel supercomputer show that the generated code meets or beats implementations that are hand-coded by a human expert and outperforms the widely used ScaLAPACK library.

Codesign Tradeoffs for High-Performance, Low-Power Linear Algebra Architectures

As technology is reaching physical limits, reducing power consumption is a key issue on our path to sustained performance. In this paper, we study fundamental tradeoffs and limits in efficiency (as measured in energy per operation) that can be achieved for an important class of kernels, namely the level-3 Basic Linear Algebra Subprograms (BLAS). It is well-accepted that specialization is the key to efficiency. This paper establishes a baseline by studying GEneral Matrix-matrix Multiplication (GEMM) on a variety of custom and general-purpose CPU and GPU architectures. Our analysis shows that orders of magnitude improvements in efficiency are possible with relatively simple customizations and fine-tuning of memory hierarchy configurations. We argue that these customizations can be generalized to perform other representative linear algebra operations. In addition to exposing the sources of inefficiencies in current CPUs and GPUs, our results show our prototype Linear Algebra Processor (LAP) implementing Double-precision GEMM (DGEMM) can achieve 600 GFLOPS while consuming less than 25 Watts in standard 45 nm technology, which is up to 50 × more energy efficient than cutting-edge CPUs.

Modeling performance through memory-stalls

We aim at modeling the performance of linear algebra algorithms without executing either them or parts of them. The performance of an algorithm can be expressed in terms of the time spent on CPU execution and on memory-stalls. The main concern of this paper is to build analytical models to accurately predict memory-stalls. The scenario in which data resides in the L2 cache is considered; with this assumption, only L1 cache misses occur. We construct an analytical formula for modeling the L1 cache misses of fundamental linear algebra operations such as those included in the Basic Linear Algebra Subprograms (BLAS) library. The number of cache misses occurring in higher-level algorithms "like a matrix factorization" is then predicted by combining the models for the appropriate BLAS subroutines. As case studies, we consider GER, a BLAS level-2 operation, and the LU factorization. The models are validated on both Intel and AMD processors, attaining remarkably accurate performance predictions.

Advanced complex trait analysis.

The Genome-wide Complex Trait Analysis (GCTA) software package can quantify the contribution of genetic variation to phenotypic variation for complex traits. However, as those datasets of interest continue to increase in size, GCTA becomes increasingly computationally prohibitive. We present an adapted version, Advanced Complex Trait Analysis (ACTA), demonstrating dramatically improved performance. We restructure the genetic relationship matrix (GRM) estimation phase of the code and introduce the highly optimized parallel Basic Linear Algebra Subprograms (BLAS) library combined with manual parallelization and optimization. We introduce the Linear Algebra PACKage (LAPACK) library into the restricted maximum likelihood (REML) analysis stage. For a test case with 8999 individuals and 279,435 single nucleotide polymorphisms (SNPs), we reduce the total runtime, using a compute node with two multi-core Intel Nehalem CPUs, from ∼17 h to ∼11 min. The source code is fully available under the GNU Public License, along with Linux binaries. For more information see http://www.epcc.ed.ac.uk/software-products/acta. a.gray@ed.ac.uk Supplementary data are available at Bioinformatics online.

A New Vectorization Technique for Expression Templates in C++

Vector operations play an important role in high performance computing and are typically provided by highly optimized libraries that implement the Basic Linear Algebra Subprograms (BLAS) interface. In C++ templates and operator overloading allow the implementation of these vector operations as expression templates which construct custom loops at compile time and providing a more abstract interface. Unfortunately existing expression template libraries lack the performance of fast BLAS implementations. This paper presents a new approach - Statically Accelerated Loop Templates (SALT) - to close this performance gap by combining expression templates with an aggressive loop unrolling technique. Benchmarks were conducted using the Intel C++ compiler and GNU Compiler Collection to assess the performance of our library relative to Intel's Math Kernel Library as well as the Eigen template library. The results show that the approach is able to provide optimization comparable to the fastest available BLAS implementations, while retaining the convenience and flexibility of a template library.

Geocomputation and Open Source Software: Components and Software Stacks

Geocomputation, with its necessary focus on software development and methods innovation, has enjoyed a close relationship with free and open source software communities. These extend from communities providing the numerical infrastructure for computation, such as BLAS (Basic Linear Algebra Subprograms), through language communities around Python, Java and others, to communities supporting spatial data handling, especially the projects of the Open Source Geospatial Foundation. This chapter surveys the stack of software components available for geocomputation from these sources, looking in most detail at the R language and environment, and how OSGeo projects have been interfaced with it. In addition, attention will be paid to open development models and community participation in software development. Since free and open source geospatial software has also achieved a successively greater presence in proprietary software as computational platforms evolve, the chapter will close with some indications of future trends in software component stacks, using Terralib as an example.

Taxonomic analysis of marine phytoplankton

Samples of sea water contain phytoplankton taxa in varying amounts, and marine scientists are interested in the relative abundance of each taxa. Their relative biomass can be ascertained indirectly by measuring the quantity of various pigments using high performance liquid chromatography. However, the conversion from pigment to taxa is mathematically non trivial as it is a positive matrix factorisation problem where both matrices are unknown beyond the level of initial estimates. The prior information on the pigment to taxa conversion matrix is used to give the problem a unique solution. An iteration of two non-negative least squares algorithms gives satisfactory results. Some sample analysis of data indicates prospects for this type of analysis. An alternative more computationally intensive approach using Bayesian methods is discussed. References Efron, B., and Tibirani, R. J., 1993. An introduction to the Bootstrap, Chapman and Hall, New York. Hastings, W. K., 1970. Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika, 57:97--109. doi:10.1093/biomet/57.1.97 Lawson C. L., and Hanson R. L. 1995. Solving Least Squares Problems, SIAM, Philadelphia, Ch 23. Lee D. D., and Seung H. S. 2001. Algorithms for Non-negative Matrix Factorization, Advances in Neural Information Processing Systems, 13:556--562. http://luthuli.cs.uiuc.edu/ daf/courses/Optimization/Papers/lee01algorithms.pdf Mackey, M. D., Mackey, D. J., Higgins, H. W. and Wright, S. W., 1996. CHEMTAX---A program for estimating class abundances from chemical markers: application to HPLC measurement of phytoplankton, Marine Ecology - Progress Series, 144:266--283. doi:10.3354/meps144265 van den Meersche, K., Soetaert, K., and Middleburg, J. J., 2008. A Bayesian computational estimator for microbial taxonomy based on biomarkers, Limnol. Oceanogr. Methods 6: 190--199. doi:10.4319/lom.2008.6.190 Robert, C. P. and Casella, G., 2004. Monte Carlo statistical methods, Springer, New York. Roberts, G. O. and Rosenthal, J. S., 2001. Optimal Scaling for Various Metropolis--Hastings Algorithms, Statistical Science 16:351--367. doi:10.1214/ss/1015346320 Wikipedia, 2011. Basic linear algebra subprograms, http://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms Wright, S. W., van den Enden, R. L., Pearce, I., Davidson, A. T. Scott, F. J., and Westwood, K. J., 2010. Phytoplankton community structure and stocks in the Southern Ocean (30-80E) determined by CHEMTAX analysis of HPLC pigment signatures Deep-Sea Research II 57:758--778. doi:10.1016/j.dsr2.2009.06.015

Speeding up solving of differential matrix Riccati equations using GPGPU computing and MATLAB

SUMMARYIn this work, we developed a parallel algorithm to speed up the resolution of differential matrix Riccati equations using a backward differentiation formula algorithm based on a fixed‐point method. The role and use of differential matrix Riccati equations is especially important in several applications such as optimal control, filtering, and estimation. In some cases, the problem could be large, and it is interesting to speed it up as much as possible. Recently, modern graphic processing units (GPUs) have been used as a way to improve performance. In this paper, we used an approach based on general‐purpose computing on graphics processing units. We used NVIDIA © GPUs with unified architecture. To do this, a special version of basic linear algebra subprograms for GPUs, called CUBLAS, and a package (three different packages were studied) to solve linear systems using GPUs have been used. Moreover, we developed a MATLAB © toolkit to use our implementation from MATLAB in such a way that if the user has a graphic card, the performance of the implementation is improved. If the user does not have such a card, the algorithm can also be run using the machine CPU. Experimental results on a NVIDIA Quadro FX 5800 are shown. Copyright © 2011 John Wiley & Sons, Ltd.

Acceleration of Hessenberg Reduction for Nonsymmetric Eigenvalue Problems in a Hybrid CPU-GPU Computing Environment

Solution of large-scale dense nonsymmetric eigenvalue problem is required in many areas of scientific and engineering computing, such as vibration analysis of automobiles and analysis of electronic diffraction patterns. In this study, we focus on the Hessenberg reduction step and consider accelerating it in a hybrid CPU-GPU computing environment. Considering that the Hessenberg reduction algorithm consists almost entirely of BLAS (Basic Linear Algebra Subprograms) operations, we propose three approaches for distributing the BLAS operations between CPU and GPU. Among them, the third approach, which assigns small-size BLAS operations to CPU and distributes large-size BLAS operations between CPU and GPU in some optimal manner, was found to be consistently faster than the other two approaches. On a machine with an Intel Core i7 processor and an NVIDIA Tesla C1060 GPU, this approach achieved 3.2 times speedup over the CPU-only case when computing the Hessenberg form of a 8,192 × 8,192 real matrix.

Scheduling Two-Sided Transformations Using Tile Algorithms on Multicore Architectures

The objective of this paper is to describe, in the context of multicore architectures, three different scheduler implementations for the two-sided linear algebra transformations, in particular the Hessenberg and Bidiagonal reductions which are the first steps for the standard eigenvalue problems and the singular value decompositions respectively. State-of-the-art dense linear algebra softwares, such as the LAPACK and ScaLAPACK libraries, suffer performance losses on multicore processors due to their inability to fully exploit thread-level parallelism. At the same time the fine-grain dataflow model gains popularity as a paradigm for programming multicore architectures. Buttari et al. (Parellel Comput. Syst. Appl. 35 (2009), 38–53) introduced the concept oftile algorithmsin which parallelism is no longer hidden inside Basic Linear Algebra Subprograms but is brought to the fore to yield much better performance. Along with efficient scheduling mechanisms for data-driven execution, these tile two-sided reductions achieve high performance computing by reaching up to 75% of the DGEMM peak on a 12000×12000 matrix with 16 Intel Tigerton 2.4 GHz processors. The main drawback of thetile algorithmsapproach for two-sided transformations is that the full reduction cannot be obtained in one stage. Other methods have to be considered to further reduce the band matrices to the required forms.

Out-of-core solution of linear systems on graphics processors

We combine two high-level application programming interfaces to solve large-scale linear systems with the data stored on disk using current graphics processors (GPU). The result is a simple yet powerful tool that enables a fast development of object-oriented codes, implemented as Matlab M-scripts, for linear algebra operations. The approach enhances the programmability of the solutions in this problem domain while unleashing the high performance of GPU. Experimental results are reported from Octave, linked with the implementation of basic linear algebra subprograms by Nvidia; as an additional novelty, we evaluate the performance of our approach for the Nvidia new double-precision graphics hardware.

High performance BLAS formulation of the adaptive Fast Multipole Method

In a previous work, we have presented a new formulation of the uniform version of the Fast Multipole Method (FMM) for the Laplace equation by using matrix products that can be efficiently computed thanks to the BLAS (Basic Linear Algebra Subprograms) routines. We propose here to extend this formulation to the adaptive version of the FMM: this requires the conception of a new data structure for the octree, namely the octree with indirections, which is efficient for both uniform and non-uniform distributions, as well as a detection mechanism of the available uniform areas in non-uniform distributions. In comparison with other M2L computation schemes (block FFT, rotations and plane wave expansions) in the case of non-uniform distributions of particles, our BLAS version appears to be the fastest for the common precisions used when one solves the Laplace equation.

C++ OPPS, a new software for the interpretation of protein dynamics from nuclear magnetic resonance measurements

AbstractNuclear magnetic resonance (NMR) is a powerful tool for elucidating protein dynamics because of the possibility to interpret nuclear spin relaxation properties in terms of microdynamic parameters. Magnetic relaxation times T1, T2, and NOE depend on dipolar and quadrupolar interactions, on chemical shift anisotropy and cross‐correlation effects. Within the framework of given motional model, it is possible to express the NMR relaxation times as functions of spectral densities (Abragam, The Principles of Nuclear Magnetism; Oxford University Press: Clarendon, London, 1961), obtaining the connection between macroscopic observables and microscopic properties. In this context, recently Meirovitch et al. (Shapiro et al., Biochemistry 2002, 41, 6271, Meirovitch et al., J Phys Chem B 2006, 110, 20615, Meirovitch et al., J Phys Chem B 2007, 111, 12865) applied the dynamical model introduced by Polimeno and Freed (Polimeno and Freed, Adv Chem Phys 1993, 83, 89, Polimeno and Freed, J Phys Chem 1995, 99, 10995), known as the slowly relaxing local structure (SRLS) model, to the study of NMR data.The program C++OPPS (http://www.chimica.unipd.it/licc/), developed in our laboratory, implements the SRLS model in an user‐friendly way with a graphical user interface (GUI), introduced to simplify the work to users who do not feel at ease with the complex mathematics of the model and the difficulties of command line based programs. The program is an evolution of the old FORTRAN 77 implementation COPPS (COupled Protein Probe Smoluchowski) and presents a number of new features: the presence of an easy to use GUI written in JAVA; high calculation performance thanks to features of C++ language, employment of BLAS (basic linear algebra subprograms) library (Blackford et al., Trans Math Soft 2002, 28, 135) in handling matrix‐vector operations and parallelization of the code under the MPI (message passing interface) paradigm (Gropp et al., Parallel Comput 1996, 22, 789, Gropp and Lusk, User's Guide for mpich, a Portable Implementation of MPI Mathematics and Computer Science Division; Argonne National Laboratory, 1996); possibility to predict the diffusion tensor of the protein via a hydrodynamic approach (Barone et al., J Comp Chem, in press). A cluster version of C++OPPS was also developed, which can be easily accessed by users via the web. © 2009 Wiley Periodicals, Inc. Int J Quantum Chem, 2010

Performance Evaluation of Multi-Core Intel Xeon Processors on Basic Linear Algebra Subprograms

Multi-core technology is a natural next step in delivering the benefits of Moore's law to computing platforms. On multi-core processors, the performance of many applications would be improved by parallel processing threads of codes using multi-threading techniques. This paper evaluates the performance of the multi-core Intel Xeon processors on the widely used basic linear algebra subprograms (BLAS). On two dual-core Intel Xeon processors with Hyper-Threading technology, our results show that a performance of around 20 GFLOPS is achieved on Level-3 (matrix-matrix operations) BLAS using multi-threading, SIMD, matrix blocking, and loop unrolling techniques. However, on a small size of Level-2 (matrix-vector operations) and Level-1 (vector operations) BLAS, the use of multi-threading technique speeds down the execution because of the thread creation overheads. Thus the use of Intel SIMD instruction set is the way to improve the performance of single-threaded Level-2 (6 GFLOPS) and Level-1 BLAS (3 GFLOPS). When the problem size becomes large (cannot fit in L2 cache), the performance of the four Xeon cores is less than 2 and 1 GFLOPS on Level-2 and Level-1 BLAS, respectively, even though eight threads are executed in parallel on eight logical processors.

HeteroPBLAS: A Set of Parallel Basic Linear Algebra Subprograms Optimized for Heterogeneous Computational Clusters

This paper presents a software library, called Heterogeneous PBLAS (HeteroPBLAS), which provides optimized parallel basic linear algebra subprograms for Heterogeneous Computational Clusters. This library is written on the top of HeteroMPI and PBLAS whose building blocks, the de facto standard kernels for matrix and vector operations (BLAS) and message passing communication (BLACS), are optimized for heterogeneous computational clusters. This is the first step towards the development of a parallel linear algebra package for Heterogeneous Computational Clusters. We show that the efficiency of the parallel routines is due to the most important feature of the library, which is the automation of the difficult optimization tasks of parallel programming on heterogeneous computing clusters. They are the determination of the accurate values of the platform parameters such as the speeds of the processors and the latencies and bandwidths of the communication links connecting different pairs of processors, the optimal values of the algorithmic parameters such as the data distribution blocking factor, the total number of processes, the 2D process grid arrangement, and the efficient mapping of the processes executing the parallel algorithm to the executing nodes of the heterogeneous computing cluster. We present the user interface and the software hierarchy of the first research implementation of HeteroPBLAS. We demonstrate the efficiency of the HeteroPBLAS programs on a homogeneous computing cluster and a heterogeneous computing cluster.