Linear Algebra Libraries Research Articles

Fast and accurate normalization of vectors and quaternions

We present fast and accurate ways to normalize two and three-dimensional vectors and quaternions and compute their length. Our approach is an adaptation of ideas used in the linear algebra library LAPACK, and we believe that the computational geometry and computer-aided design communities are not aware of the possibility of speeding up these fundamental operations in the robust way proposed here.

Incompressible SPH (ISPH) with fast Poisson solver on a GPU

This paper presents a fast incompressible SPH (ISPH) solver implemented to run entirely on a graphics processing unit (GPU) capable of simulating several millions of particles in three dimensions on a single GPU. The ISPH algorithm is implemented by converting the highly optimised open-source weakly-compressible SPH (WCSPH) code DualSPHysics to run ISPH on the GPU, combining it with the open-source linear algebra library ViennaCL for fast solutions of the pressure Poisson equation (PPE). Several challenges are addressed with this research: constructing a PPE matrix every timestep on the GPU for moving particles, optimising the limited GPU memory, and exploiting fast matrix solvers. The ISPH pressure projection algorithm is implemented as 4 separate stages, each with a particle sweep, including an algorithm for the population of the PPE matrix suitable for the GPU, and mixed precision storage methods. An accurate and robust ISPH boundary condition ideal for parallel processing is also established by adapting an existing WCSPH boundary condition for ISPH. A variety of validation cases are presented: an impulsively started plate, incompressible flow around a moving square in a box, and dambreaks (2-D and 3-D) which demonstrate the accuracy, flexibility, and speed of the methodology. Fragmentation of the free surface is shown to influence the performance of matrix preconditioners and therefore the PPE matrix solution time. The Jacobi preconditioner demonstrates robustness and reliability in the presence of fragmented flows. For a dambreak simulation, GPU speed ups demonstrate up to 10–18 times and 1.1–4.5 times compared to single-threaded and 16-threaded CPU run times respectively.

Базовые библиотеки линейной алгебры для высокопроизводительных расчетов

The article considers linear algebra libraries such as BLAS, LAPACK, ScaLAPACK, MKL, and ATLAS, which support high-performance computing (HPC) in modern architectures and are used both in well-known performance tests and in various applications. In the majority of applications, the most time-consuming computation stages are implemented by calling subroutines from such libraries; therefore, the optimal choice of a library is an important issue in setting up computations. The main aim of this review is to describe the invariant characteristics of libraries to achieve high performance of applications. High performance computations used in different fields of knowledge are briefly reviewed. Classification of linear algebra libraries in terms of their functionality and applied high-performance architectures is suggested. The basic low-level BLAS library implemented for all HPC architectures is demonstrated. It is pointed out that the BLAS library supports dividing of the entire computation process into several parallel flows in systems with a common memory field; for such systems, tools such as OpenMP or OpenACC are used. In the case of systems with distributed memory, the parallel version of this library, called PBLAS is used, which supports exchange of messages between nodes using the MPI standard. Higher-level libraries based on the BLAS, e.g., the LAPACK library, which contains a large set of different programs for linear algebra, are described. The ScaLAPACK library for the distributed memory model, which is based on the LAPACK and PBLAS libraries, as well as the Intel MKL library, which is its modern development, are presented. To support efficient operation of hybrid systems, the fundamentally new libraries MAGMA and PLASMA involving features for optimizing linear-algebraic computing of small dimension are analyzed. Libraries supporting solution of eigenvalue problems, such as the EISPACK, PeigS, and a number of other libraries, are investigated. It is pointed out that in the new ELPA library oriented to supercomputers, both OpenMP and MPI tools can be used. It is noted that operations on sparse matrices, especially multiplication of matrices, are very relevant for many applied fields of science; in this regard, the SparseBLAS library can be considered to be the basic standard for them. It is concluded that the optimal choice of a library depends essentially on both the particular application and on the used computing architecture.

Using Spherical-Harmonics Expansions for Optics Surface Reconstruction from Gradients.

In this paper, we propose a new algorithm to reconstruct optics surfaces (aka wavefronts) from gradients, defined on a circular domain, by means of the Spherical Harmonics. The experimental results indicate that this algorithm renders the same accuracy, compared to the reconstruction based on classical Zernike polynomials, using a smaller number of polynomial terms, which potentially speeds up the wavefront reconstruction. Additionally, we provide an open-source C++ library, released under the terms of the GNU General Public License version 2 (GPLv2), wherein several polynomial sets are coded. Therefore, this library constitutes a robust software alternative for wavefront reconstruction in a high energy laser field, optical surface reconstruction, and, more generally, in surface reconstruction from gradients. The library is a candidate for being integrated in control systems for optical devices, or similarly to be used in ad hoc simulations. Moreover, it has been developed with flexibility in mind, and, as such, the implementation includes the following features: (i) a mock-up generator of various incident wavefronts, intended to simulate the wavefronts commonly encountered in the field of high-energy lasers production; (ii) runtime selection of the library in charge of performing the algebraic computations; (iii) a profiling mechanism to measure and compare the performance of different steps of the algorithms and/or third-party linear algebra libraries. Finally, the library can be easily extended to include additional dependencies, such as porting the algebraic operations to specific architectures, in order to exploit hardware acceleration features.

Algorithm 979

To exploit both memory locality and the full performance potential of highly tuned kernels, dense linear algebra libraries, such as linear algebra package (LAPACK), commonly implement operations as blocked algorithms. However, to achieve near-optimal performance with such algorithms, significant tuning is required. In contrast, recursive algorithms are virtually tuning free and attain similar performance. In this article, we first analyze and compare blocked and recursive algorithms in terms of performance and then introduce recursive LAPACK (R e LAPACK), an open-source library of recursive algorithms to seamlessly replace many of LAPACK’s blocked algorithms. In most scenarios, R e LAPACK outperforms reference LAPACK and in many situations improves upon the performance of optimized libraries.

Time and space efficient generators for quasiseparable matrices

The class of quasiseparable matrices is defined by the property that any submatrix entirely below or above the main diagonal has small rank, namely below a bound called the order of quasiseparability. These matrices arise naturally in solving PDE's for particle interaction with the Fast Multi-pole Method (FMM), or computing generalized eigenvalues. From these application fields, structured representations and algorithms have been designed in numerical linear algebra to compute with these matrices in time linear in the matrix dimension and either quadratic or cubic in the quasiseparability order. Motivated by the design of the general purpose exact linear algebra library LinBox, and by algorithmic applications in algebraic computing, we adapt existing techniques introduce novel ones to use quasiseparable matrices in exact linear algebra, where sub-cubic matrix arithmetic is available. In particular, we will show, the connection between the notion of quasiseparability and the rank profile matrix invariant, that we have introduced in 2015. It results in two new structured representations, one being a simpler variation on the hierarchically semiseparable storage, and the second one exploiting the generalized Bruhat decomposition. As a consequence, most basic operations, such as computing the quasiseparability orders, applying a vector, a block vector, multiplying two quasiseparable matrices together, inverting a quasiseparable matrix, can be at least as fast and often faster than previous existing algorithms.

A high-performance Riccati based solver for tree-structured quadratic programs

Robust multi-stage Model Predictive Control (MPC) is an increasingly popular approach to handle model uncertainties due to the simplicity of its problem formulation and other attractive properties. However, the exponential growth of the problem dimensions with respect to the robust horizon renders the online solution of such problems challenging and the development of tailored solvers crucial. In this paper, an interior point method is presented that can solve Quadratic Programs (QPs) arising in multi-stage MPC efficiently by means of a tree-structured Riccati recursion and a high-performance linear algebra library. A performance comparison with code-generated and general purpose sparse QP solvers shows that the computation times can be significantly reduced for all problem sizes that are practically relevant in embedded MPC applications. The presented implementation is freely available as part of the open-source software HPMPC.

SMArtOp: A Java library for distributing high-dimensional sparse-matrix arithmetic operations

• Sparse-matrix operations are commonplace in computational science. • Existing libraries are mostly unable to cope with high-dimensional data. • SMArtOp is a Java library for distributing large-scale sparse-matrix arithmetic operations. • SMArtOp relies on the characteristics of operations and the matrices involved. • SMArtOp improved execution times and resource needs of popular algebra libraries. Sparse-matrix operations are commonplace in computational science, and novel solutions for speeding-up them are essential for numerous applications. SMArtOp is a software for efficiently dividing and distributing the processing of large-scale sparse-matrix arithmetic operations. This software relies on both the distinctive characteristics of each type of arithmetic operation and the particular matrices involved to split the operations into parallel and simpler tasks. Experimental evaluation showed the speeding-up and resource consumption advantages of the proposed software, in comparison to other linear-algebra libraries.

Cholesky factorization on SIMD multi-core architectures

Abstract Many linear algebra libraries, such as the Intel MKL, Magma or Eigen, provide fast Cholesky factorization. These libraries are suited for big matrices but perform slowly on small ones. Even though State-of-the-Art studies begin to take an interest in small matrices, they usually feature a few hundreds rows. Fields like Computer Vision or High Energy Physics use tiny matrices. In this paper we show that it is possible to speed up the Cholesky factorization for tiny matrices by grouping them in batches and using highly specialized code. We provide High Level Transformations that accelerate the factorization for current multi-core and many-core SIMD architectures ( SSE, AVX 2, KNC, AVX512 , Neon, Altivec). We focus on the fact that, on some architectures, compilers are unable to vectorize and on other architectures, vectorizing compilers are not efficient. Thus hand-made SIMD ization is mandatory. We achieve with these transformations combined with SIMD a speedup from × 14 to × 28 for the whole resolution in single precision compared to the naive code on a AVX 2 machine and a speedup from × 6 to × 14 on double precision, both with a strong scalability.

Programming With Models: Writing Statistical Algorithms for General Model Structures With NIMBLE

ABSTRACTWe describe NIMBLE, a system for programming statistical algorithms for general model structures within R. NIMBLE is designed to meet three challenges: flexible model specification, a language for programming algorithms that can use different models, and a balance between high-level programmability and execution efficiency. For model specification, NIMBLE extends the BUGS language and creates model objects, which can manipulate variables, calculate log probability values, generate simulations, and query the relationships among variables. For algorithm programming, NIMBLE provides functions that operate with model objects using two stages of evaluation. The first stage allows specialization of a function to a particular model and/or nodes, such as creating a Metropolis-Hastings sampler for a particular block of nodes. The second stage allows repeated execution of computations using the results of the first stage. To achieve efficient second-stage computation, NIMBLE compiles models and functions via C++, using the Eigen library for linear algebra, and provides the user with an interface to compiled objects. The NIMBLE language represents a compilable domain-specific language (DSL) embedded within R. This article provides an overview of the design and rationale for NIMBLE along with illustrative examples including importance sampling, Markov chain Monte Carlo (MCMC) and Monte Carlo expectation maximization (MCEM). Supplementary materials for this article are available online.

Tuning linear algebra for energy efficiency on multicore machines by adapting the ATLAS library

While automated tuning is an established method for minimising the execution time of scientific applications, it has rarely been used for an automated minimisation of the energy consumption. This article presents a study on how to adapt the auto-tuned linear algebra library ATLAS to consider the energy consumption of the execution in its tuning decision. For different tuning parameters of ATLAS, it investigates which differences occur in the tuning results when ATLAS is tuned for a minimal execution time or for a minimal energy consumption. The tuning parameters include the matrix size for the low-level matrix multiplication, loop unrolling factors, crossover points for different matrix-multiplication implementations, the minimum size for matrices to be transposed, or blocking sizes for the last-level cache. Also, parameters for multithreaded execution, such as the number of threads and thread affinity are investigated. The emphasis of this article is on a method proposed with which it is possible to replace a tuning process for execution time by a tuning for energy consumption, especially in the parallel case. ATLAS serves as a prominent example for a tuned library. Furthermore, the article draws conclusions on how to design an energy-optimising autotuning package and how to choose tuning parameters. The article also discusses why the matrix-matrix multiplication has a potential for increasing the energy efficiency while the time efficiency remains constant, whereas other routines have shown to improve their energy efficiency by reducing the execution time.

Task-based Cholesky decomposition on Xeon Phi architectures using OpenMP

The increasing number of computational cores in modern many-core processors, as represented by the Intel Xeon Phi architectures, has created the need for an open-source, high performance and scalable task-based dense linear algebra package that can efficiently use this type of many-core hardware. In this paper, we examined the design modifications necessary when porting PLASMA, a task-based dense linear algebra library, run effectively on two generations of Intel's Xeon Phi architecture, known as knights corner (KNC) and knights landing (KNL). First, we modified PLASMA's tiled Cholesky decomposition to use OpenMP tasks for its scheduling mechanism to enable Xeon Phi compatibility. We then compared the performance of our modified code to that of the original dynamic scheduler running on an Intel Xeon Sandy Bridge CPU. Finally, we looked at the performance of the OpenMP tiled Cholesky decomposition on knights corner and knights landing processors. We detail the optimisations required to obtain performance on these platforms and compare with the highly tuned Intel MKL math library.

Multi-threaded dense linear algebra libraries for low-power asymmetric multicore processors

Dense linear algebra libraries, such as BLAS and LAPACK, provide a relevant collection of numerical tools for many scientific and engineering applications. While there exist high performance implementations of the BLAS (and LAPACK) functionality for many current multi-threaded architectures, the adaption of these libraries for asymmetric multicore processors (AMPs) is still pending. In this paper we address this challenge by developing an asymmetry-aware implementation of the BLAS, based on the BLIS framework, and tailored for AMPs equipped with two types of cores: fast/power-hungry versus slow/energy-efficient. For this purpose, we integrate coarse-grain and fine-grain parallelization strategies into the library routines which, respectively, dynamically distribute the workload between the two core types and statically repartition this work among the cores of the same type.Our results on an ARM® big.LITTLE™ processor embedded in the Exynos 5422 SoC, using the asymmetry-aware version of the BLAS and a plain migration of the legacy version of LAPACK, experimentally assess the benefits, limitations, and potential of this approach from the perspectives of both throughput and energy efficiency.

An implementation of the direct-forcing immersed boundary method using GPU power

ABSTRACTA graphics processing unit (GPU) is utilized to apply the direct-forcing immersed boundary method. The code running on the GPU is generated with the help of the Compute Unified Device Architecture (CUDA). The first and second spatial derivatives of the incompressible Navier-Stokes equations are discretized by the sixth-order central compact finite-difference schemes. Two flow fields are simulated. The first test case is the simulated flow around a square cylinder, with the results providing good estimations of the wake region mechanics and vortex shedding. The second test case is the simulated flow around a circular cylinder. This case was selected to better understand the effects of sharp corners on the force coefficients. It was observed that the estimation of the force coefficients did not result in any troubles in the case of a circular cylinder. Additionally, the performance values obtained for the calculation time for the solution of the Poisson equation are compared with the values for other CPUs and GPUs from the literature. Consequently, approximately 3× and 20× speedups are achieved in comparison with GPU (using CUSP library) and CPU, respectively. CUSP is an open-source library for sparse linear algebra and graph computations on CUDA.

GPU-based parallel method of temperature field analysis in a floor heater with a controller

Abstract A parallel method enabling acceleration of the numerical analysis of the transient temperature field in an air floor heating system is presented in this paper. An initial-boundary value problem of the heater regulated by an on/off controller is formulated. The analogue model is discretized using the implicit finite difference method. The BiCGStab method is used to compute the obtained system of equations. A computer program implementing simultaneous computations on CPUand GPU(GPGPUtechnology) was developed. CUDA environment and linear algebra libraries (CUBLAS and CUSPARSE) are used by this program. The time of computations was reduced eight times in comparison with a program executed on the CPU only. Results of computations are presented in the form of time profiles and temperature field distributions. An influence of a model of the heat transfer coefficient on the simulation of the system operation was examined. The physical interpretation of obtained results is also presented.Results of computations were verified by comparing them with solutions obtained with the use of a commercial program - COMSOL Mutiphysics.

Armadillo: a template-based C++ library for linear algebra

The C++ language is often used for implementing functionality that is performance and/or resource sensitive. While the standard C++ library provides many useful algorithms (such as sorting), in its current form it does not provide direct handling of linear algebra (matrix maths). Armadillo is an open source linear algebra library for the C++ language, aiming towards a good balance between speed and ease of use. Its high-level Application Programming Interface (API) is deliberately similar to the widely Matlab and Octave languages (Eaton et al. 2015), so that mathematical operations can be expressed in a familiar and natural manner. The library is useful for algorithm development directly in C++, or relatively quick conversion of research code into production environments. Armadillo provides efficient objects for vectors, matrices and cubes (third order tensors), as well as over 200 associated functions for manipulating data stored in the objects. Integer, floating point and complex numbers are supported, as well as dense and sparse storage formats. Various matrix factorisations are provided through integration with LAPACK (Demmel 1997), or one of its high performance drop-in replacements such as Intel MKL (Intel 2016) or OpenBLAS (Xianyi, Qian, and Saar 2016). It is also possible to use Armadillo in conjunction with NVBLAS to obtain GPU-accelerated matrix multiplication (NVIDIA 2015). Armadillo is used as a base for other open source projects, such as MLPACK, a C++ library for machine learning and pattern recognition (Curtin et al. 2013), and RcppArmadillo, a bridge between the R language and C++ in order to speed up computations (Eddelbuettel and Sanderson 2014). Armadillo internally employs an expression evaluator based on template metaprogramming techniques (Abrahams and Gurtovoy 2004), to automatically combine several operations in order to increase speed and efficiency. An overview of the internal architecture is given in (Sanderson 2010).

The BLIS Framework

BLIS is a new software framework for instantiating high-performance BLAS-like dense linear algebra libraries. We demonstrate how BLIS acts as a productivity multiplier by using it to implement the level-3 BLAS on a variety of current architectures. The systems for which we demonstrate the framework include state-of-the-art general-purpose, low-power, and many-core architectures. We show, with very little effort, how the BLIS framework yields sequential and parallel implementations that are competitive with the performance of ATLAS, OpenBLAS (an effort to maintain and extend the GotoBLAS), and commercial vendor implementations such as AMD’s ACML, IBM’s ESSL, and Intel’s MKL libraries. Although most of this article focuses on single-core implementation, we also provide compelling results that suggest the framework’s leverage extends to the multithreaded domain.

LinearOperator—A generic, high-level expression syntax for linear algebra

We introduce an expression syntax for the evaluation of matrix–matrix, matrix–vector and vector–vector operations. The implementation is similar to the well-known general concept of expression templates as used, for example, in the C++ linear-algebra libraries Eigen and Blaze. The novelty of the approach that is discussed here lies in the use of new C++11 features like lambda expressions and std::function objects that avoid the majority of the implementational complexity that usually comes with a pure template solution.A concrete implementation of the expression syntax has been developed within the framework of the finite-element library deal.II, but it is fairly generic: the LinearOperator implementation only requires a minimal vector and matrix interface, that all of deal.II’s concrete vector and matrix types adhere to. This makes the interface fully transparent with respect to the concrete implementation, in particular to the storage strategy (full matrix, sparse structure), and memory strategy (local, shared, distributed).The paper concludes with a number of performance comparisons and examples that demonstrate that the framework results in efficient, short and concise code. The performance comparisons show that the overhead introduced by std::function objects is negligible for moderately sized matrices, even when compared to native expression-template implementations.

Parallel Programming Models for Dense Linear Algebra on Heterogeneous Systems

We present a review of the current best practices in parallel programming models for dense linear algebra DLA on heterogeneous architectures. We consider multicore CPUs, stand alone manycore coprocessors, GPUs, and combinations of these. Of interest is the evolution of the programming models for DLA libraries { in particular, the evolution from the popular LAPACK and ScaLAPACK libraries to their modernized counterparts PLASMA for multicore CPUs and MAGMA for heterogeneous architectures, as well as other programming models and libraries. Besides providing insights into the programming techniques of the libraries considered, we outline our view of the current strengths and weaknesses of their programming models { especially in regards to hardware trends and ease of programming high-performance numerical software that current applications need { in order to motivate work and future directions for the next generation of parallel programming models for high-performance linear algebra libraries on heterogeneous systems.

PyVCI: A flexible open-source code for calculating accurate molecular infrared spectra

The PyVCI program package is a general purpose open-source code for simulating accurate molecular spectra, based upon force field expansions of the potential energy surface in normal mode coordinates. It includes harmonic normal coordinate analysis and vibrational configuration interaction (VCI) algorithms, implemented primarily in Python for accessibility but with time-consuming routines written in C. Coriolis coupling terms may be optionally included in the vibrational Hamiltonian. Non-negligible VCI matrix elements are stored in sparse matrix format to alleviate the diagonalization problem. CPU and memory requirements may be further controlled by algorithmic choices and/or numerical screening procedures, and recommended values are established by benchmarking using a test set of 44 molecules for which accurate analytical potential energy surfaces are available.Force fields in normal mode coordinates are obtained from the PyPES library of high quality analytical potential energy surfaces (to 6th order) or by numerical differentiation of analytic second derivatives generated using the GAMESS quantum chemical program package (to 4th order). Program summaryProgram title: PyVCICatalogue identifier: AFAC_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AFAC_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: MIT LicenseNo. of lines in distributed program, including test data, etc.: 401703No. of bytes in distributed program, including test data, etc.: 7091446Distribution format: tar.gzProgramming language: Python, C.Computer: PC.Operating system: Linux, MacOSX, Windows.RAM: Varies widelyClassification: 16.3.External routines: Numpy, Scipy, CythonNature of problem: The simulation of accurate molecular vibrational spectra is a significant and long-standing problem in computational chemistry. There are two major challenges: constructing an accurate ab initio potential energy surface and solving the nuclear vibrational problem. Both scale poorly with respect to molecular size, requiring large amounts of CPU time and memory.Solution method: We have implemented a straightforward numerical differentiation algorithm to construct quartic force fields in normal mode coordinates using second derivatives of the energy with respect to nuclear displacement obtained from ab initio quantum chemical calculations, for nuclear vibrational structure algorithm development and testing purposes. We have also provided an interface to the PyPES library of high quality semi-global potential energy surfaces, which enable quantitative prediction of molecular vibrational spectra. To solve the nuclear vibrational problem, we use a vibrational configuration interaction algorithm in a harmonic oscillator basis.Unusual features: One of the unusual features of our code is its flexibility, with multiple ways of generating or supplying force field data, dynamic memory allocation, adjustable screening thresholds, and explicit user control over terms in the VCI wave-function (maximum excitation level and extent of mode-coupling). We employ sparse matrix linear algebra libraries to reduce the memory required for VCI matrix storage and diagonalization, and provide for parallel VCI matrix construction to reduce required wall times.Additional comments: User Manual and examples (tutes) includedRunning time: Varies widely