Linear algebra subprograms on shared memory computers

The modern graphics processing unit (GPU) found in many standard personal computers is a highly parallel math processor capable of nearly 1 TFLOPS peak throughput at a cost similar to a high-end CPU and an excellent FLOPS/watt ratio. High-level linear algebra operations are computationally intense, often requiring O(N3) operations and would seem a natural fit for the processing power of the GPU. Our work is on CULA, a GPU accelerated implementation of linear algebra routines. We present results from factorizations such as LU decomposition, singular value decomposition and QR decomposition along with applications like system solution and least squares. The GPU execution model featured by NVIDIA GPUs based on CUDA demands very strong parallelism, requiring between hundreds and thousands of simultaneous operations to achieve high performance. Some constructs from linear algebra map extremely well to the GPU and others map poorly. CPUs, on the other hand, do well at smaller order parallelism and perform acceptably during low-parallelism code segments. Our work addresses this via hybrid a processing model, in which the CPU and GPU work simultaneously to produce results. In many cases, this is accomplished by allowing each platform to do the work it performs most naturally.

The aim of this communication is to present a Symbolic Linear Algebra package, written in MAPLE V.5 and working also in the release 7, with general but especially educational purpose. Its goal is to run over MAPLE the different Linear Algebra algorithms developped in [OJ]. The implemented algorithms are also valids over finite characteristic and, for educational purpose, may be runned interactively. Three topics are covered by the package: equivalent matrices (echelon matrices, rank, LU decomposition, linear systems,...), similar matrices (Eigenvalues and eigenvectors, Rational or Frobenius form, Irreducible form, Jordan form,...), and congruent matrices (Symmetric and Hermitian matrices, Cholesky and QR decomposition, SVD, Gram-Schmidt over arbitrary euclidean metrics, ortogonal matrices, quadric surfaces...)

This paper is an updated and extended version of the paper "The QR Decomposition and the Singular Value Decomposition in the Symmetrized Max-Plus Algebra" (B. De Schutter and B. De Moor, SIAM J. Matrix Anal. Appl., 19 (1998), pp. 378--406). The max-plus algebra, which has maximization and addition as its basic operations, can be used to describe and analyze certain classes of discrete-event systems, such as flexible manufacturing systems, railway networks, and parallel processor systems. In contrast to conventional algebra and conventional (linear) system theory, the max-plus algebra and the max-plus-algebraic system theory for discrete-event systems are at present far from fully developed, and many fundamental problems still have to be solved. Currently, much research is going on to deal with these problems and to further extend the max-plus algebra and to develop a complete max-plus-algebraic system theory for discrete-event systems. In this paper we address one of the remaining gaps in the max-plus algebra by considering matrix decompositions in the symmetrized max-plus algebra. The symmetrized max-plus algebra is an extension of the max-plus algebra obtained by introducing a max-plus-algebraic analogue of the $-$-operator. We show that we can use well-known linear algebra algorithms to prove the existence of max-plus-algebraic analogues of basic matrix decompositions from linear algebra such as the QR decomposition, the singular value decomposition, the Hessenberg decomposition, the LU decomposition, and so on. These max-plus-algebraic matrix decompositions could play an important role in the max-plus-algebraic system theory for discrete-event systems.

We present algorithms to determine the number of nonzeros in each row and column of the factors of a sparse matrix, for both the QR factorization and the LU factorization with partial pivoting. The algorithms use only the nonzero structure of the input matrix, and run in time nearly linear in the number of nonzeros in that matrix. They may be used to set up data structures or schedule parallel operations in advance of the numerical factorization. The row and column counts we compute are upper bounds on the actual counts. If the input matrix is strong Hall and there is no coincidental numerical cancellation, the counts are exact for QR factorization and are the tightest bounds possible for LU factorization. These algorithms are based on our earlier work on computing row and column counts for sparse Cholesky factorization, plus an efficient method to compute the column elimination tree of a sparse matrix without explicitly forming the product of the matrix and its transpose.

QR and LU decompositions are the most important matrix decomposition algorithms. Many studies work on accelerating these algorithms by FPGA or ASIC in a case by case style. In this paper, we propose a unified framework for the matrix decomposition algorithms, combining three QR decomposition algorithms and LU algorithm with pivoting into a unified linear array structure. The QR and LU decomposition algorithms exhibit the same two-level loop structure and the same data dependency. Utilizing the similarities in loop structure and data dependency of matrix decomposition, we unify a fine-grained algorithm for all four matrix decomposition algorithms. Furthermore, we present a unified coprocessor structure with a scalable linear array of processing elements (PEs), in which four types of PEs are same in the structure of memory channels and PE connections, but the only difference exists in the internal structure of data path. Our unified co-processor, which is IEEE 32-bit floating-point precision, is implemented and mapped onto a Xilinx Virtex5 FPGA chip. Experimental results show that our co-processors can achieve speedup of 2.3 to 14.9 factors compared to a Pentium Dual CPU with double SSE threads.

Chapter 9 - Numerical Techniques

1. Linear systems: Introduction to Linear Systems. Gauss Elimination. Numerical Solutions. 2. Vectors: Vector Operations. Dot Product. Span. Linear Independence. The Product [Ax]. Cross Product. Lines, Planes, and Hyperplanes. 3. Matrices: Matrix Operations. Matrix Inverse. Elementary and Invertible Matrices. LU Factorization. 4. Vector spaces: Subspaces of [R to the n Power]. Vector Spaces. Linear Independence Bases. Dimension. Coordinate Vectors and Change of Basis. Rank and Nullity. Applications to Coding Theory. 5. Linear transformations: Matrix Transformations. Linear Transformations. Kernel and Range. The Matrix of a Linear Transformation. The Algebra of Linear Transformations. 6. Determinants: Determinants Cofactor Expansion. Properties of Determinants. The Adjoint Cramer's Rule. Determinants with Permutations. 7. Eigenvalues and eigenvectors: Eigenvalues and Eigenvectors. Diagonalization. Approximations of Eigenvalues and Eigenvectors. Applications to Dynamical Systems. Applications to Markov Chains. 8. Dot and inner products: Orthogonal Sets and Matrices. Orthogonal Projections Gram-Schmidt Process. The QR Factorization. Least Squares. Orthogonalization of Symmetric Matrices. Quadratic Forms and Conic Sections. The Singular Value Decomposition (SVD). Inner Products. Appendices: Answers to selected exercises Each chapter concludes with Applications, Mini-Projects, and Computer Exercises.

This paper intends to show the impact of compiler optimization flags and the variable's precision on direct methods to solve linear systems. The chosen six methods are simple direct methods, so our work could be a study for new researchers in this field. The methods are LU decomposition, LDU decomposition, Gaussian Elimination, Gauss-Jordan Elimination, Cholesky decomposition, and QR decomposition using the Gram-Schmidt orthogonalization process. Our study showed a huge difference in time between single- and double-precision in all methods, but the error encountered in the single-precision was not so high. Also, the best flags to these methods were the `-O3' and the `-Ofast' ones.

In Gaussian elimination and LU factorization we solve a linear system by transforming it to triangular form. These are not the only kind of transformations that can be used for such a task. Matrices with orthonormal columns, called unitary matrices can be used to reduce a square matrix to upper triangular form and more generally a rectangular matrix to upper triangular (also called upper trapezoidal) form. This lead to a decomposition of a rectangular matrix known as a QR decomposition and a reduced form which we refer to as a QR factorization. The QR decomposition and factorization will be used in later chapters to solve least squares- and eigenvalue problems.

From the Publisher: This is the revised and greatly expanded Second Edition of the hugely popular Numerical Recipes: The Art of Scientific Computing. The product of a unique collaboration among four leading scientists in academic research and industry, Numerical Recipes is a complete text and reference book on scientific computing. In a self-contained manner it proceeds from mathematical and theoretical considerations to actual practical computer routines. With over 100 new routines (now well over 300 in all), plus upgraded versions of many of the original routines, this book is more than ever the most practical, comprehensive handbook of scientific computing available today. The book retains the informal, easy-to-read style that made the first edition so popular, with many new topics presented at the same accessible level. In addition, some sections of more advanced material have been introduced, set off in small type from the main body of the text. Numerical Recipes is an ideal textbook for scientists and engineers and an indispensable reference for anyone who works in scientific computing. Highlights of the new material include a new chapter on integral equations and inverse methods; multigrid methods for solving partial differential equations; improved random number routines; wavelet transforms; the statistical bootstrap method; a new chapter on less-numerical algorithms including compression coding and arbitrary precision arithmetic; band diagonal linear systems; linear algebra on sparse matrices; Cholesky and QR decomposition; calculation of numerical derivatives; Pade approximants, and rational Chebyshev approximation; new special functions; Monte Carlo integration in high-dimensional spaces; globally convergent methods for sets of nonlinear equations; an expanded chapter on fast Fourier methods; spectral analysis on unevenly sampled data; Savitzky-Golay smoothing filters; and two-dimensional Kolmogorov-Smirnoff tests. All this is in addition to material on such basic top

A class of parallel tiled linear algebra algorithms for multicore architectures

The sparse matrix solver has become the bottleneck in a Simulation Program with Integrated Circuit Emphasis circuit simulator. It is difficult to parallelize the sparse matrix solver because of the high data dependence during the numerical LU factorization. In this brief, a parallel LU factorization algorithm is developed on shared-memory computers with multicore central processing units, based on KLU algorithms. An Elimination Scheduler (EScheduler) is proposed to represent the data dependence during the LU factorization. Based on the EScheduler, the parallel tasks are scheduled in two modes to achieve a high level of concurrence, i.e., cluster mode and pipeline mode . The experimental results on 26 circuit matrices reveal that the developed algorithm can achieve speedup of 1.18-4.55t (on geometric average), as compared with KLU, with 1-8 threads. The result analysis indicates that for different data dependence, different parallel strategies should be dynamically selected to obtain optimal performance.

Dimension reduction plays an important role in handling the massive quantity of high dimensional data such as biomedical text data, gene expression data, and mass spectrometry data, and so forth. In this paper, we introduce distance preserving dimension reduction (DPDR) based on the QR factorization (DPDR/QR) or the Cholesky factorization (DPDR/C). DPDR generates lower dimensional representations of the high-dimensional data, which can exactly preserve Euclidean distances and cosine similarities between any pair of data points in the original dimensional space. After projecting data points to the lower dimensional space obtained from DPDR, one can execute other data analysis algorithms. DPDR can substantially reduce the computing time and/or memory requirement of a given data analysis algorithm, especially when we need to run the data analysis algorithm many times for estimating parameters or searching for a better solution.

The article investigates the use of linear algebra methods for query optimisation in modern database management systems. The relevance of the study is due to the growth of data volumes and the need for more efficient methods of their processing. The proposed approach is to formalise tabular data structures as numerical matrices, which allows us to interpret relational operations (selection, projection, join) as linear algebraic transformations, which opens up opportunities for the use of powerful mathematical tools that optimise computations. The main methods of linear algebra for query optimisation are considered: matrix factorisation (LU, QR decomposition) for structuring complex calculations; singular value decomposition (SVD) for reducing dimensionality and identifying hidden dependencies; matrix modelling of JOIN operations as an efficient block multiplication; systems of linear equations for formalising data filtering. It is noted that matrix factorisation helps to optimise the execution of complex analytical queries, SVD is effective in aggregation and dimensionality reduction tasks, the matrix representation of JOIN provides high performance when processing large amounts of data, and linear systems unify the filtering logic with numerical processing. To empirically verify the proposed methods, an experimental study was conducted on a synthetic dataset simulating a transactional table with a million records. The experiment was carried out using the NumPy and PyTorch libraries with the use of GPU acceleration. We evaluated such indicators as query execution time, RAM usage, and scalability. The results showed that JOIN matrix modelling was the fastest, which is especially important for parallel computing. LU/QR factorisation demonstrated stability and versatility, while SVD proved to be resource-intensive but valuable for data compression. Linear equation methods showed balanced performance. It is concluded that linear algebra is a powerful formalisation framework for next-generation DBMS architectures capable of providing mathematical consistency, computational efficiency and scalability.

In Demmel et al. (Numer. Math. 106(2), 199–224, 2007) we showed that a large class of fast recursive matrix multiplication algorithms is stable in a normwise sense, and that in fact if multiplication of n-by-n matrices can be done by any algorithm in O(n ω+η ) operations for any η > 0, then it can be done stably in O(n ω+η ) operations for any η > 0. Here we extend this result to show that essentially all standard linear algebra operations, including LU decomposition, QR decomposition, linear equation solving, matrix inversion, solving least squares problems, (generalized) eigenvalue problems and the singular value decomposition can also be done stably (in a normwise sense) in O(n ω+η ) operations.