Linear Algebra Routines Research Articles

X-OpenMP — eXtreme fine-grained tasking using lock-less work stealing

Processors with 100s of threads of execution are among the state-of-the-art in high-end computing systems. This transition to many-core computing has required the community to develop new algorithms to overcome significant latency bottlenecks through massive concurrency. However, implementing efficient parallel runtimes that can scale up to high concurrency levels with extremely fine-grained tasks remains a challenge. Existing techniques do not scale to a large number of threads due to the high cost of synchronization in concurrent data structures. We present a thorough analysis of various synchronization mechanisms including mutex, semaphore, spinlock and atomic fetch-and-add that are typically used to build concurrent data structures in task-parallel runtime systems. To overcome these limitations, in a recent work we proposed XQueue, a novel lock-less concurrent queuing system with relaxed ordering semantics that is geared towards realizing scalability up to hundreds of concurrent threads. In this work, we extend XQueue and present X-OpenMP, a library for enabling extremely fine-grained parallelism on modern many-core systems with hundreds of cores. Work stealing is a popular choice for load balancing in task-based runtime systems as it efficiently distributes the load across worker threads; however, traditional approaches rely on synchronization primitives and thus work stealing can incur overheads. Here we implement a lock-less algorithm for work stealing for total-store order (TSO) memory architectures and evaluate the performance using micro and macro benchmarks. We compare the performance of X-OpenMP with native LLVM OpenMP, GNU OpenMP, OpenCilk and oneTBB implementations using task-based linear algebra routines from PLASMA numerical library, Strassen’s matrix multiplication from the BOTS Benchmark Suite, and the Unbalanced Tree Search benchmark. Applications parallelized using OpenMP can run without modification by simply linking against the X-OpenMP library. X-OpenMP achieves up to 40X speedup compared to GNU OpenMP, up to 2X speedup compared to the native LLVM OpenMP, up to 6X speedup compared to OpenCilk and up to 5X speedup compared to oneTBB implementations. The tasking overheads in X-OpenMP are reduced by 50% compared to the native LLVM OpenMP.

BLAS Kütüphanelerinin GPU Mimarilerindeki Nicel Performans Analizi

Basic Linear Algebra Subprograms (BLAS) are a set of linear algebra routines commonly used by machine learning applications and scientific computing. BLAS libraries with optimized implementations of BLAS routines offer high performance by exploiting parallel execution units in target computing systems. With massively large number of cores, graphics processing units (GPUs) exhibit high performance for computationally-heavy workloads. Recent BLAS libraries utilize parallel cores of GPU architectures efficiently by employing inherent data parallelism. In this study, we analyze GPU-targeted functions from two BLAS libraries, cuBLAS and MAGMA, and evaluate their performance on a single-GPU NVIDIA architecture by considering architectural features and limitations. We collect architectural performance metrics and explore resource utilization characteristics. Our work aims to help researchers and programmers to understand the performance behavior and GPU resource utilization of the BLAS routines implemented by the libraries.

Efficient computation of optical excitations in two-dimensional materials with the Xatu code

Here we describe an efficient numerical implementation of the Bethe-Salpeter equation to obtain the excitonic spectrum of semiconductors. This is done on the electronic structure calculated either at the simplest tight-binding level or through density functional theory calculations based on local orbitals. We use a simplified model for the electron-electron interactions which considers atomic orbitals as point-like orbitals and a phenomenological screening. The optical conductivity can then be optionally computed within the Kubo formalism. Our results for paradigmatic two-dimensional materials such as hBN and MoS2, when compared with those of more sophisticated first-principles methods, are excellent and envision a practical use of our implementation beyond the computational limitations of such methods. Program summaryProgram Title: XatuCPC Library link to program files:https://doi.org/10.17632/kj4rt95pvc.1Developer's repository link:https://github.com/alejandrojuria/xatuLicensing provisions: GPLv3Programming language: C++, Fortran, PythonNature of problem: The exciton spectrum is obtained as the solution of the Bethe-Salpeter equation for insulators and semi-conductors. Constructing the equation involves determining the screening of the electrostatic interaction and then determining the matrix elements of the interaction kernel, which are computationally-intensive tasks, specially if one takes a purely ab-initio approach.Solution method: The Bethe-Salpeter equation can be efficiently set up and solved assuming that the basis of the reference electronic structure calculation, obtained either from tight-binding models or density functional theory with actual localized orbitals, corresponds to point-like localized orbitals. This, in addition to using an effective screening instead of computing the dielectric constant, allows to obtain the interaction kernel at very low computational cost and, thereof, the exciton spectrum as well as the light absorption of materials.Additional comments including restrictions and unusual features: The code requires using at least C++11, given that it uses version-specific features. All linear algebra routines have been delegated to the Armadillo library.

A quadratic decoder approach to nonintrusive reduced‐order modeling of nonlinear dynamical systems

AbstractLinear projection schemes like Proper Orthogonal Decomposition can efficiently reduce the dimensions of dynamical systems but are naturally limited, e.g., for convection‐dominated problems. Nonlinear approaches have shown to outperform linear methods in terms of dimension reduction versus accuracy but, typically, come with a large computational overhead. In this work, we consider a quadratic reduction scheme which induces nonlinear structures that are well accessible to tensorized linear algebra routines. We discuss that nonintrusive approaches can be used to simultaneously reduce the complexity in the equations and propose an operator inference formulation that respects dynamics on nonlinear manifolds.

Proximal Stabilized Interior Point Methods and Low-Frequency-Update Preconditioning Techniques

In this work, in the context of Linear and convex Quadratic Programming, we consider Primal Dual Regularized Interior Point Methods (PDR-IPMs) in the framework of the Proximal Point Method. The resulting Proximal Stabilized IPM (PS-IPM) is strongly supported by theoretical results concerning convergence and the rate of convergence, and can handle degenerate problems. Moreover, in the second part of this work, we analyse the interactions between the regularization parameters and the computational footprint of the linear algebra routines used to solve the Newton linear systems. In particular, when these systems are solved using an iterative Krylov method, we are able to show—using a new rearrangement of the Schur complement which exploits regularization—that general purposes preconditioners remain attractive for a series of subsequent IPM iterations. Indeed, if on the one hand a series of theoretical results underpin the fact that the approach here presented allows a better re-use of such computed preconditioners, on the other, we show experimentally that such (re)computations are needed only in a fraction of the total IPM iterations. The resulting regularized second order methods, for which low-frequency-update of the preconditioners are allowed, pave the path for an alternative class of second order methods characterized by reduced computational effort.

The triple decomposition of the velocity gradient tensor as a standardized real Schur form

The triple decomposition of a velocity gradient tensor provides an analysis tool in fluid mechanics by which the flow can be split into a sum of irrotational straining flow, shear flow, and rigid body rotational flow. In 2007, Kolář formulated an optimization problem to compute the triple decomposition [V. Kolář, “Vortex identification: New requirements and limitations,” Int. J. Heat Fluid Flow 28, 638–652 (2007)], and more recently, the triple decomposition has been connected to the Schur form of the associated matrix. We show that the standardized real Schur form, which can be computed by state of the art linear algebra routines, is a solution to the optimization problem posed by Kolář. We also demonstrate why using the standardized variant of the real Schur form makes computation of the triple decomposition more efficient. Furthermore, we illustrate why different structures of the real Schur form correspond to different alignments of the coordinate system with the fluid flow and may, therefore, lead to differences in the resulting triple decomposition. Based on these results, we propose a new, simplified algorithm for computing the triple decomposition, which guarantees consistent results.

Providing performance portable numerics for Intel GPUs

SummaryWith discrete Intel GPUs entering the high‐performance computing landscape, there is an urgent need for production‐ready software stacks for these platforms. In this article, we report how we enable the Ginkgo math library to execute on Intel GPUs by developing a kernel backed based on the DPC++ programming environment. We discuss conceptual differences between the CUDA and DPC++ programming models and describe workflows for simplified code conversion. We evaluate the performance of basic and advanced sparse linear algebra routines available in Ginkgo's DPC++ backend in the hardware‐specific performance bounds and compare against routines providing the same functionality that ship with Intel's oneMKL vendor library.

Discrete Lehmann representation of imaginary time Green's functions

We present an efficient basis for imaginary time Green's functions based on a low-rank decomposition of the spectral Lehmann representation. The basis functions are simply a set of well-chosen exponentials, so the corresponding expansion may be thought of as a discrete form of the Lehmann representation using an effective spectral density which is a sum of $\ensuremath{\delta}$ functions. The basis is determined only by an upper bound on the product $\ensuremath{\beta}{\ensuremath{\omega}}_{max}$, with $\ensuremath{\beta}$ the inverse temperature and ${\ensuremath{\omega}}_{max}$ an energy cutoff, and a user-defined error tolerance $\ensuremath{\epsilon}$. The number $r$ of basis functions scales as $O(log(\ensuremath{\beta}{\ensuremath{\omega}}_{max})log(1/\ensuremath{\epsilon}))$. The discrete Lehmann representation of a particular imaginary time Green's function can be recovered by interpolation at a set of $r$ imaginary time nodes. Both the basis functions and the interpolation nodes can be obtained rapidly using standard numerical linear algebra routines. Due to the simple form of the basis, the discrete Lehmann representation of a Green's function can be explicitly transformed to the Matsubara frequency domain, or obtained directly by interpolation on a Matsubara frequency grid. We benchmark the efficiency of the representation on simple cases, and with a high-precision solution of the Sachdev-Ye-Kitaev equation at low temperature. We compare our approach with the related intermediate representation method, and introduce an improved algorithm to build the intermediate representation basis and a corresponding sampling grid.

QPALM: a proximal augmented lagrangian method for nonconvex quadratic programs

We propose QPALM, a nonconvex quadratic programming (QP) solver based on the proximal augmented Lagrangian method. This method solves a sequence of inner subproblems which can be enforced to be strongly convex and which therefore admit a unique solution. The resulting steps are shown to be equivalent to inexact proximal point iterations on the extended-real-valued cost function, which allows for a fairly simple analysis where convergence to a stationary point at an \(R\)-linear rate is shown. The QPALM algorithm solves the subproblems iteratively using semismooth Newton directions and an exact linesearch. The former can be computed efficiently in most iterations by making use of suitable factorization update routines, while the latter requires the zero of a monotone, one-dimensional, piecewise affine function. QPALM is implemented in open-source C code, with tailored linear algebra routines for the factorization in a self-written package LADEL. The resulting implementation is shown to be extremely robust in numerical simulations, solving all of the Maros-Meszaros problems and finding a stationary point for most of the nonconvex QPs in the Cutest test set. Furthermore, it is shown to be competitive against state-of-the-art convex QP solvers in typical QPs arising from application domains such as portfolio optimization and model predictive control. As such, QPALM strikes a unique balance between solving both easy and hard problems efficiently.

A fast data-driven method for genotype imputation, phasing and local ancestry inference: MendelImpute.jl.

Current methods for genotype imputation and phasing exploit the volume of data in haplotype reference panels and rely on hidden Markov models (HMMs). Existing programs all have essentially the same imputation accuracy, are computationally intensive and generally require prephasing the typed markers. We introduce a novel data-mining method for genotype imputation and phasing that substitutes highly efficient linear algebra routines for HMM calculations. This strategy, embodied in our Julia program MendelImpute.jl, avoids explicit assumptions about recombination and population structure while delivering similar prediction accuracy, better memory usage and an order of magnitude or better run-times compared to the fastest competing method. MendelImpute operates on both dosage data and unphased genotype data and simultaneously imputes missing genotypes and phase at both the typed and untyped SNPs (single nucleotide polymorphisms). Finally, MendelImpute naturally extends to global and local ancestry estimation and lends itself to new strategies for data compression and hence faster data transport and sharing. Software, documentation and scripts to reproduce our results are available from https://github.com/OpenMendel/MendelImpute.jl. Supplementary data are available at Bioinformatics online.

A Heuristic Independent Particle Approximation to Determinantal Point Processes

A determinantal point process is a stochastic point process that is commonly used to capture negative correlations. It has become increasingly popular in machine learning in recent years. Sampling a determinantal point process however remains a computationally intensive task. This note introduces a heuristic independent particle approximation to determinantal point processes. The approximation is based on the physical intuition of fermions and is implemented using standard numerical linear algebra routines. Sampling from this independent particle approximation can be performed at a negligible cost. Numerical results are provided to demonstrate the performance of the proposed algorithm.

A survey of numerical linear algebra methods utilizing mixed-precision arithmetic

The efficient utilization of mixed-precision numerical linear algebra algorithms can offer attractive acceleration to scientific computing applications. Especially with the hardware integration of low-precision special-function units designed for machine learning applications, the traditional numerical algorithms community urgently needs to reconsider the floating point formats used in the distinct operations to efficiently leverage the available compute power. In this work, we provide a comprehensive survey of mixed-precision numerical linear algebra routines, including the underlying concepts, theoretical background, and experimental results for both dense and sparse linear algebra problems.

Design and implementation of a modular interior-point solver for linear optimization

This paper introduces the algorithmic design and implementation of Tulip, an open-source interior-point solver for linear optimization. It implements a regularized homogeneous interior-point algorithm with multiple centrality corrections, and therefore handles unbounded and infeasible problems. The solver is written in Julia, thus allowing for a flexible and efficient implementation: Tulip’s algorithmic framework is fully disentangled from linear algebra implementations and from a model’s arithmetic. In particular, this allows to seamlessly integrate specialized routines for structured problems. Extensive computational results are reported. We find that Tulip is competitive with open-source interior-point solvers on the H. Mittelmann’s benchmark of barrier linear programming solvers. Furthermore, we design specialized linear algebra routines for structured master problems in the context of Dantzig–Wolfe decomposition. These routines yield a tenfold speedup on large and dense instances that arise in power systems operation and two-stage stochastic programming, thereby outperforming state-of-the-art commercial interior point method solvers. Finally, we illustrate Tulip’s ability to use different levels of arithmetic precision by solving problems in extended precision.

Efficient Auto-Tuning of Parallel Programs with Interdependent Tuning Parameters via Auto-Tuning Framework (ATF)

Auto-tuning is a popular approach to program optimization: it automatically finds good configurations of a program’s so-called tuning parameters whose values are crucial for achieving high performance for a particular parallel architecture and characteristics of input/output data. We present three new contributions of the Auto-Tuning Framework (ATF), which enable a key advantage in general-purpose auto-tuning : efficiently optimizing programs whose tuning parameters have interdependencies among them. We make the following contributions to the three main phases of general-purpose auto-tuning: (1) ATF generates the search space of interdependent tuning parameters with high performance by efficiently exploiting parameter constraints; (2) ATF stores such search spaces efficiently in memory, based on a novel chain-of-trees search space structure; (3) ATF explores these search spaces faster, by employing a multi-dimensional search strategy on its chain-of-trees search space representation. Our experiments demonstrate that, compared to the state-of-the-art, general-purpose auto-tuning frameworks, ATF substantially improves generating, storing, and exploring the search space of interdependent tuning parameters, thereby enabling an efficient overall auto-tuning process for important applications from popular domains, including stencil computations, linear algebra routines, quantum chemistry computations, and data mining algorithms.

Study of Exploiting Coarse‐Grained Parallelism in Block‐Oriented Numerical Linear Algebra Routines

AbstractWe have developed streaming implementations of two numerical linear algebra operations that further exploit the block decomposition strategies commonly used in these operations to obtain performance. The implementations formulate algorithms as data flow graphs and use coarse‐grained parallelism to (1) emit a block in the result matrix as soon as it becomes available and (2) compute on multiple blocks in parallel. This streaming design benefits data flow graphs consisting of multiple linear algebra operations as it removes synchronization points between successive operations: a result block from an operation can be used immediately in an algorithm's successor operations without waiting for the full result from the first operation. Early comparisons with OpenBLAS functions on CPUs show comparable performance for computing with large dense matrices and an earliest arrival time of a result block that is up to 50x smaller than the time needed for a full result. More thorough studies can show the impact of such implementations on the performance of systems by chaining multiple linear algebra operations.

Integrating software and hardware hierarchies in an autotuning method for parallel routines in heterogeneous clusters

A hierarchical approach for autotuning linear algebra routines on heterogeneous platforms is presented. Hierarchy helps to alleviate the difficulties of tuning parallel routines for high-performance computing systems. This paper analyzes the application of the hierarchical approach at both the hardware and software levels, using the basic matrix multiplication and the Strassen multiplication as proof of concept on multicore+coprocessor nodes. In this way, the hierarchical approach allows partial delegation of the efficient exploitation of the computing units in the node to the underlying direct autotuned matrix multiplication used in the base case.

Thermal-FIST: A package for heavy-ion collisions and hadronic equation of state

Thermal-FIST11Thermal-FIST — Thermal, Fast and Interactive Statistical Toolkit. is a C++ package designed for convenient general-purpose physics analysis within the family of hadron resonance gas (HRG) models. This mainly includes the statistical analysis of particle production in heavy-ion collisions and the phenomenology of hadronic equation of state. Notable features include fluctuations and correlations of conserved charges, effects of probabilistic decay, chemical non-equilibrium, and inclusion of van der Waals hadronic interactions. Calculations are possible within the grand canonical ensemble, the canonical ensemble, as well as in mixed-canonical ensembles combining the canonical treatment of certain conserved charges with the grand-canonical treatment of other conserved charges. The package contains a fast thermal event generator, which generates particle yields in accordance with the HRG chemistry, and particle momenta based on the Blast Wave model. A distinct feature of this package is the presence of the graphical user interface frontend – QtThermalFIST – which is designed for fast and convenient general-purpose HRG model applications. Program summaryProgram Title:Thermal-FIST, version 1.2Program Files doi:http://dx.doi.org/10.17632/pprr8p4fkp.1Licensing provisions: GPLv3Programming language: C++External routines:Eigen template library for the linear algebra routines [1], MINUIT2 package from CERN ROOT [2], Mersenne Twister random number generator [3], Qt5 framework [4] (for the GUI only), QCustomPlot Qt widget [5] (for the GUI only)Nature of problem: The HRG model and its various modifications constitute a common framework used for modeling of the hadronic equation of state and particle production in heavy-ion collisions. Even the simplest versions of the HRG model require careful considerations of the many details, including the resonance decay feed-down, implementation of charge conservation constraints relevant for heavy-ion collisions, chemical non-equilibrium effects. A notable extra effort is required in order to treat the fluctuations and correlations of various charges , which presently are being extensively studied in the heavy-ion collision experiments and lattice QCD calculations. The inclusion of hadronic interactions, modeled by an excluded-volume (EV) or a van der Waals (vdW) type framework, additionally requires a numerical solution to a system of many transcendental equations.Solution method: The Thermal-FIST package contains a class-based library which calculates relevant HRG observables for a specified setup. The setup includes a particle list, usually to be supplied with an external file, an HRG model specification (statistical ensemble, van der Waals interaction parameters, etc.), a set of thermal parameters, and conservation laws constraints. Whenever necessary, the systems of transcendental equations are solved numerically with the Broyden’s method. The package includes a fitter for extracting thermal parameters from hadron yield data through the χ2 minimization. The HRG model based Monte Carlo event generator is a complementary feature to analytic calculations. A general-purpose thermal analysis is made maximally convenient with QtThermalFIST — a GUI frontend based on the Qt framework where all typical calculations, such as the properties of the equation of state or the thermal fits, can be straightforwardly performed.Additional comments: If the EV/vdW interactions are present, exact analytic calculations are presently only possible within the grand canonical ensemble. Approximate calculations are possible for the strangeness-canonical ensemble on the condition that strange particles form a small subsystem relative to the total system. Effects of probabilistic decays on fluctuation observables are generally included only up to the moments of the 2nd order. The only exception is the ideal HRG model in the grand canonical ensemble, where these effects are included up to the moments of the 4th order. On the other hand, the Monte Carlo event generator is not constrained by the above restrictions.

BLASFEO

Basic Linear Algebra Subroutines for Embedded Optimization (BLASFEO) is a dense linear algebra library providing high-performance implementations of BLAS- and LAPACK-like routines for use in embedded optimization and small-scale high-performance computing, in general. A key difference with respect to existing high-performance implementations of BLAS is that the computational performance is optimized for small- to medium-scale matrices, i.e., for sizes up to a few hundred. BLASFEO comes with three different implementations: a high-performance implementation aimed at providing the highest performance for matrices fitting in cache, a reference implementation providing portability and embeddability and optimized for very small matrices, and a wrapper to standard BLAS and LAPACK providing high performance on large matrices. The three implementations of BLASFEO together provide high-performance dense linear algebra routines for matrices ranging from very small to large. Compared to both open-source and proprietary highly tuned BLAS libraries, for matrices of size up to about 100, the high-performance implementation of BLASFEO is about 20--30% faster than the corresponding level 3 BLAS routines and two to three times faster than the corresponding LAPACK routines.

Modal Analysis of Fluid Flows: An Overview

Modal Analysis of Fluid Flows: An Overview

Fast synchronization‐free algorithms for parallel sparse triangular solves with multiple right‐hand sides

SummaryThe sparse triangular solve kernels, SpTRSV and SpTRSM, are important building blocks for a number of numerical linear algebra routines. Parallelizing SpTRSV and SpTRSM on today's manycore platforms, such as GPUs, is not an easy task since computing a component of the solution may depend on previously computed components, enforcing a degree of sequential processing. As a consequence, most existing work introduces a preprocessing stage to partition the components into a group of level‐sets or colour‐sets so that components within a set are independent and can be processed simultaneously during the subsequent solution stage. However, this class of methods requires a long preprocessing time as well as significant runtime synchronization overheads between the sets. To address this, we propose in this paper novel approaches for SpTRSV and SpTRSM in which the ordering between components is naturally enforced within the solution stage. In this way, the cost for preprocessing can be greatly reduced, and the synchronizations between sets are completely eliminated. To further exploit the data‐parallelism, we also develop an adaptive scheme for efficiently processing multiple right‐hand sides in SpTRSM. A comparison with a state‐of‐the‐art library supplied by the GPU vendor, using 20 sparse matrices on the latest GPU device, shows that the proposed approach obtains an average speedup of over two for SpTRSV and up to an order of magnitude speedup for SpTRSM. In addition, our method is up to two orders of magnitude faster for the preprocessing stage than existing SpTRSV and SpTRSM methods.