Dense Matrix Operations Research Articles

Distributed out-of-memory NMF on CPU/GPU architectures

AbstractWe propose an efficient distributed out-of-memory implementation of the non-negative matrix factorization (NMF) algorithm for heterogeneous high-performance-computing systems. The proposed implementation is based on prior work on NMFk, which can perform automatic model selection and extract latent variables and patterns from data. In this work, we extend NMFk by adding support for dense and sparse matrix operation on multi-node, multi-GPU systems. The resulting algorithm is optimized for out-of-memory problems where the memory required to factorize a given matrix is greater than the available GPU memory. Memory complexity is reduced by batching/tiling strategies, and sparse and dense matrix operations are significantly accelerated with GPU cores (or tensor cores when available). Input/output latency associated with batch copies between host and device is hidden using CUDA streams to overlap data transfers and compute asynchronously, and latency associated with collective communications (both intra-node and inter-node) is reduced using optimized NVIDIA Collective Communication Library (NCCL) based communicators. Benchmark results show significant improvement, from 32X to 76x speedup, with the new implementation using GPUs over the CPU-based NMFk. Good weak scaling was demonstrated on up to 4096 multi-GPU cluster nodes with approximately 25,000 GPUs when decomposing a dense 340 Terabyte-size matrix and an 11 Exabyte-size sparse matrix of density $$10^{-6}$$ 10 - 6 .

Unleashing the Potential of Sparse DNNs Through Synergistic Hardware-Sparsity Co-Design

Sparsity, a widely recognized path to curbing the computational needs of deep neural networks (DNNs), still suffers a number of roadblocks in practice, despite a decade of intensive research on sparse neural networks. The extant structured sparsity patterns often fail to attain significant model compression, while the hardware challenges posed by unstructured sparsity are yet to be fully overcome. As algorithmic and hardware innovations individually deliver limited benefits, a synergistic approach is necessary to unleash the potential of sparse DNNs. This work proposes a tightly integrated design methodology for the sparsity patterns and associated hardware platforms to reach the highest model compression goals while simultaneously facilitating efficient hardware processing. We demonstrate that novel complementary sparsity patterns can offer utmost expressiveness levels with inherent hardware exploitable regularity. Our novel dynamic training method converts the expressiveness of such sparsity configurations into highly accurate and compact sparse neural networks. Complementary sparsity is represented in a dense format, and when synergistically coupled with minimal yet strategic hardware modifications, can be processed in close concordance with the conventional dataflow of the dense matrix operations. We thus demonstrate that there is ample room for innovation beyond conventional techniques and immense practical potential for sparse neural networks through the synergistic design of sparsity patterns and hardware architectures.

Very fast finite element Poisson solvers on lower precision accelerator hardware: A proof of concept study for Nvidia Tesla V100

Recently, accelerator hardware in the form of graphics cards including Tensor Cores, specialized for AI, has significantly gained importance in the domain of high-performance computing. For example, NVIDIA’s Tesla V100 promises a computing power of up to 125 TFLOP/s achieved by Tensor Cores, but only if half precision floating point format is used. We describe the difficulties and discrepancy between theoretical and actual computing power if one seeks to use such hardware for numerical simulations, that is, solving partial differential equations with a matrix-based finite element method, with numerical examples. If certain requirements, namely low condition numbers and many dense matrix operations, are met, the indicated high performance can be reached without an excessive loss of accuracy. A new method to solve linear systems arising from Poisson’s equation in 2D that meets these requirements, based on “prehandling” by means of hier-archical finite elements and an additional Schur complement approach, is presented and analyzed. We provide numerical results illustrating the computational performance of this method and compare it to a commonly used (geometric) multigrid solver on standard hardware. It turns out that we can exploit nearly the full computational power of Tensor Cores and achieve a significant speed-up compared to the standard methodology without losing accuracy.

Towards electronic structure-based ab-initio molecular dynamics simulations with hundreds of millions of atoms

We push the boundaries of electronic structure-based ab-initio molecular dynamics (AIMD) beyond 100 million atoms. This scale is otherwise barely reachable with classical force-field methods or novel neural network and machine learning potentials. We achieve this breakthrough by combining innovations in linear-scaling AIMD, efficient and approximate sparse linear algebra, low and mixed-precision floating-point computation on GPUs, and a compensation scheme for the errors introduced by numerical approximations. The core of our work is the non-orthogonalized local submatrix method (NOLSM), which scales very favorably to massively parallel computing systems and translates large sparse matrix operations into highly parallel, dense matrix operations that are ideally suited to hardware accelerators. We demonstrate that the NOLSM method, which is at the center point of each AIMD step, is able to achieve a sustained performance of 324 PFLOP/s in mixed FP16/FP32 precision corresponding to an efficiency of 67.7% when running on 1536 NVIDIA A100 GPUs.

Tensorox: Accelerating GPU Applications via Neural Approximation on Unused Tensor Cores

Driven by the demands of deep learning, many hardware accelerators, including GPUs, have begun to include specialized tensor processing units to accelerate matrix operations. However, general-purpose GPU applications that have little or no large dense matrix operations cannot benefit from these tensor units. This article proposes Tensorox, a framework that exploits the half-precision tensor cores available on recent GPUs for approximable, non deep learning applications. In essence, a shallow neural network is trained based on the input-output mapping of the function to be approximated. The key innovation in our implementation is the use of the small and dimension-restricted tensor operations in Nvidia GPUs to run multiple instances of the approximation neural network in parallel. With the proper scaling and training methods, our approximation yielded an overall accuracy that is higher than naïvely running the original programs with half-precision. Furthermore, Tensorox allows for the runtime adjustment of the degree of approximation. For the 10 benchmarks we tested, we achieved speedups from 2× to 112× compared to the original in single precision floating point, while maintaining the error caused by the approximation to below 10 percent in most applications.

The Second Competition on Spatial Statistics for Large Datasets

In the last few decades, the size of spatial and spatio-temporal datasets in many research areas has rapidly increased with the development of data collection technologies. As a result, classical statistical methods in spatial statistics are facing computational challenges. For example, the kriging predictor in geostatistics becomes prohibitive on traditional hardware architectures for large datasets as it requires high computing power and memory footprint when dealing with large dense matrix operations. Over the years, various approximation methods have been proposed to address such computational issues, however, the community lacks a holistic process to assess their approximation efficiency. To provide a fair assessment, in 2021, we organized the first competition on spatial statistics for large datasets, generated by our ExaGeoStat software, and asked participants to report the results of estimation and prediction. Thanks to its widely acknowledged success and at the request of many participants, we organized the second competition in 2022 focusing on predictions for more complex spatial and spatio-temporal processes, including univariate nonstationary spatial processes, univariate stationary space-time processes, and bivariate stationary spatial processes. In this paper, we describe in detail the data generation procedure and make the valuable datasets publicly available for a wider adoption. Then, we review the submitted methods from fourteen teams worldwide, analyze the competition outcomes, and assess the performance of each team.

Multicolor low‐rank preconditioner for general sparse linear systems

SummaryThis article presents a multilevel parallel preconditioning technique for solving general large sparse linear systems of equations. Subdomain coloring is invoked to reorder the coefficient matrix by multicoloring the adjacency graph of the subdomains, resulting in a two‐level block diagonal structure. A full binary tree structure is then built to facilitate the construction of the preconditioner. A key property that is exploited is the observation that the difference between the inverse of the original matrix and that of its block diagonal approximation is often well approximated by a low‐rank matrix. This property and the block diagonal structure of the reordered matrix lead to a multicolor low‐rank (MCLR) preconditioner. The construction procedure of the MCLR preconditioner follows a bottom‐up traversal of the tree . All irregular matrix computations, such as ILU factorizations and related triangular solves, are restricted to leaf nodes where these operations can be performed independently. Computations in nonleaf nodes only involve easy‐to‐optimize dense matrix operations. In order to further reduce the number of iteration of the Preconditioned Krylov subspace procedure, we combine MCLR with a few classical block‐relaxation techniques. Numerical experiments on various test problems are proposed to illustrate the robustness and efficiency of the proposed approach for solving large sparse symmetric and nonsymmetric linear systems.

GPU Acceleration of Dense Matrix and Block Operations for Lanczos Method for Systems Over GF(2)

GPU Acceleration of Dense Matrix and Block Operations for Lanczos Method for Systems Over GF(2)

Large-scale electromagnetic field analyses of coils wound with coated conductors using a current-vector-potential formulation with a thin-strip approximation

We developed a novel software for large-scale electromagnetic field analyses of coils wound with coated conductors based on current-vector-potential formulation with thin-strip approximation. Although this formulation was effective for obtaining the precise solutions of the electromagnetic field, the strong nonlinear property of superconducting materials frequently led to highly ill-conditioned linear systems of equations, which were difficult to solve efficiently. Moreover, the memory consumption and computation time required for the analyses rapidly increased with the size of the analysis due to dense matrix operations. In our software, the first difficulty was addressed by a novel preconditioning technique based on the algebraic multigrid method. Algebraic multigrid preconditioning enabled us to efficiently and stably solve the ill-conditioned linear systems of equations encountered in our analyses. It also improved the robustness of the analyses containing multifilament-coated conductors. As regards the second difficulty, the hierarchical matrices representation drastically reduced the memory consumption related to the dense matrices, as well as computation time. Meanwhile, our implementation of the hierarchical matrices representation was quite compatible with parallel computations on distributed memory computers. Finally, we presented some practical examples of large-scale analyses, which became possible by using the new software. For instance, the analysis of a cosine-theta dipole magnet whose number of degrees of freedom was more than 1.5 million was successfully completed in 78 h by 56 parallel processes and with a total memory consumption of 177 GB.

Network Coding in Heterogeneous Multicore IoT Nodes With DAG Scheduling of Parallel Matrix Block Operations

Random linear network coding (RLNC) has the potential to improve the performance of current and future Internet of Things (IoT) communication systems, but is computationally demanding due to matrix multiplications and inversions. Some single-core RLNC implementations achieve already sufficient coding speeds for contemporary multimedia streaming formats. However, advances in multimedia streaming formats and IoT applications will require the exploitation of heterogeneous multicore architectures, which are becoming common for a wide range of IoT nodes, including smartphones. In this paper, we introduce and evaluate efficient RLNC computing strategies for IoT node architectures, including the emerging heterogeneous big.LITTLE multicore architectures with multiple big (fast) cores and multiple LITTLE (slow) cores. In contrast to existing RLNC implementation strategies, we build on and adapt highly optimized dense matrix operations from the high performance computing field to RLNC on heterogeneous multicore IoT nodes. Our approach includes the optimization of RLNC matrix operations through optimized operations on matrix blocks with single instruction multiple data instructions. We schedule block operations on the heterogeneous cores through a directed acyclic graph that avoids artificial synchronization points while ensuring the data dependencies. We examine priority scheduling according to the number of outgoing dependencies of a task and data locality of cached blocks. Our extensive measurements with several heterogeneous big.LITTLE multicore IoT node and smartphone processor boards demonstrate higher RLNC encoding and decoding throughputs than existing approaches. Moreover, our measurements indicate that the utilization of more cores decreases energy consumption, which is an important goal for IoT nodes.

MARE2DEM: a 2-D inversion code for controlled-source electromagnetic and magnetotelluric data

This work presents MARE2DEM, a freely available code for 2-D anisotropic inversion of magnetotelluric (MT) data and frequency-domain controlled-source electromagnetic (CSEM) data from onshore and offshore surveys. MARE2DEM parametrizes the inverse model using a grid of arbitrarily shaped polygons, where unstructured triangular or quadrilateral grids are typically used due to their ease of construction. Unstructured grids provide significantly more geometric flexibility and parameter efficiency than the structured rectangular grids commonly used by most other inversion codes. Transmitter and receiver components located on topographic slopes can be tilted parallel to the boundary so that the simulated electromagnetic fields accurately reproduce the real survey geometry. The forward solution is implemented with a goal-oriented adaptive finite-element method that automatically generates and refines unstructured triangular element grids that conform to the inversion parameter grid, ensuring accurate responses as the model conductivity changes. This dual-grid approach is significantly more efficient than the conventional use of a single grid for both the forward and inverse meshes since the more detailed finite-element meshes required for accurate responses do not increase the memory requirements of the inverse problem. Forward solutions are computed in parallel with a highly efficient scaling by partitioning the data into smaller independent modeling tasks consisting of subsets of the input frequencies, transmitters and receivers. Non-linear inversion is carried out with a new Occam inversion approach that requires fewer forward calls. Dense matrix operations are optimized for memory and parallel scalability using the ScaLAPACK parallel library. Free parameters can be bounded using a new non-linear transformation that leaves the transformed parameters nearly the same as the original parameters within the bounds, thereby reducing non-linear smoothing effects. Data balancing normalization weights for the joint inversion of two or more data sets encourages the inversion to fit each data type equally well. A synthetic joint inversion of marine CSEM and MT data illustrates the algorithm's performance and parallel scaling on up to 480 processing cores. CSEM inversion of data from the Middle America Trench offshore Nicaragua demonstrates a real world application. The source code and MATLAB interface tools are freely available at http://mare2dem.ucsd.edu.

A sparse stiff chemistry solver based on dynamic adaptive integration for efficient combustion simulations

A sparse stiff chemistry solver based on dynamic adaptive hybrid integration (AHI-S) is developed and demonstrated for efficient combustion simulations. In a previous study, a dynamic adaptive method for hybrid integration (AHI) was developed to speed up the time integration of chemically reacting flows with detailed chemistry. The AHI method solves the fast subcomponent of chemistry implicitly and the slow subcomponent of chemistry and transport explicitly, and it was shown that AHI is more accurate and efficient than the operator-splitting schemes when there are significant radical sources from the transport term. In the present study, the AHI method is first improved to minimize the number of nontrivial entries in the Jacobian. Sparse matrix techniques are further integrated into AHI to achieve high computational efficiency. The performance of the new AHI-S solver is investigated in constant-pressure auto-ignition systems using different mechanisms that consist of 9–2878 species. It is shown that the computational cost of the AHI-S solver is overall linearly proportional to the mechanism size and is comparable to that of evaluating reaction rates using CHEMKIN-II subroutines. The AHI-S solver achieves speed-up factors ranging from approximately 10, for the 9-species hydrogen mechanism, to approximately 3000, for the 2878-species biodiesel mechanism, compared with the fully implicit VODE solver with Jacobian evaluated through numerical perturbations and factorized with dense matrix operations. It is further found that for mechanisms with less than approximately 100 species, the time saving of AHI-S is primarily attributed to the reduced size of the implicit core of the governing equations, while for mechanisms with more than 100 species, the computational cost of VODE is dominated by the dense LU factorization, such that the time saving of AHI-S is mostly attributed to the sparse LU factorization. The AHI-S solver is then applied to unsteady perfectly stirred reactors involving extinction and re-ignition. Speed-up factors from 50 to 30,000 are achieved compared with the Strang splitting scheme with the chemistry substeps implicitly integrated with VODE, while speed-up factors of 10–100 are achieved compared with the Strang splitting scheme implemented with the sparse stiff LSODES solver. In the end, the performance of AHI-S is investigated in one-dimensional (1-D) unsteady freely propagating laminar premixed flames for a methane/air mixture, for which the time step size in AHI-S is limited by the fastest transport process. A speed-up factor of approximately 200 is achieved compared with the Strang splitting scheme for fixed time step sizes between 10−8s and 10−6s.

A fast, memory efficient and robust sparse preconditioner based on a multifrontal approach with applications to finite‐element matrices

SummaryIn this article, we introduce a fast, memory efficient and robust sparse preconditioner that is based on a direct factorization scheme for sparse matrices arising from the finite‐element discretization of elliptic partial differential equations. We use a fast (but approximate) multifrontal approach as a preconditioner and use an iterative scheme to achieve a desired accuracy. This approach combines the advantages of direct and iterative schemes to arrive at a fast, robust, and accurate preconditioner. We will show that this approach is faster (∼2×) and more memory efficient (∼2–3×) than a conventional direct multifrontal approach. Furthermore, we will demonstrate that this preconditioner is both faster and more effective than other preconditioners such as the incomplete LU preconditioner. Specific speedups depend on the matrix size and improve as the size of the matrix increases. The preconditioner can be applied to both structured and unstructured meshes in a similar manner. We build on our previous work and utilize the fact that dense frontal and update matrices, in the multifrontal algorithm, can be represented as hierarchically off‐diagonal low‐rank matrices. Using this idea, we replace all large dense matrix operations in the multifrontal elimination process with O(N) hierarchically off‐diagonal low‐rank operations to arrive at a faster and more memory efficient factorization scheme. We then use this direct factorization method at low accuracies as a preconditioner and apply it to various real‐life engineering test cases. Copyright © 2016 John Wiley & Sons, Ltd.

Implementation of a generalized exponential basis functions method for linear and non‐linear problems

SummaryIn this paper, we address shortcomings of the method of exponential basis functions by extending it to general linear and non‐linear problems. In linear problems, the solution is approximated using a linear combination of exponential functions. The coefficients are calculated such that the homogenous form of equation is satisfied on some grid. To solve non‐linear problems, they are converted to into a succession of linear ones using a Newton–Kantorovich approach. The generalized exponential basis functions (GEBF) method developed can be implemented with greater ease compared with exponential basis functions, as all calculations can be performed using real numbers and no characteristic equation is needed. The details of an optimized implementation are described. We compare GEBF on some benchmark problems with methods in the literature, such as variants of the boundary element method, where GEBF shows a good performance. Also, in a 3D problem, we report the run time of the proposed method compared with that of Kratos, a parallel, highly optimized finite element code. The results show that in this example, to obtain the same level of error, much less computational effort is needed in the proposed method. Practical limitations might be encountered, however, for large problems because of dense matrix operations involved. Copyright © 2015 John Wiley & Sons, Ltd.

Faster solvers for large kinetic mechanisms using adaptive preconditioners

The adaptive preconditioners developed in this paper substantially reduce the computational cost of integrating large kinetic mechanisms using implicit ordinary differential equation (ODE) solvers. For a well-stirred reactor, the speedup of the new method is an order of magnitude faster than recent approaches based on direct, sparse linear system solvers. Moreover, the new method is up to three orders of magnitude faster than traditional implementations of the ODE solver where the Jacobian information is generated automatically via finite differences, and the factorization relies on standard, dense matrix operations. Unlike mechanism reduction strategies, the adaptive preconditioners do not alter the underlying system of differential equations. Consequently, the new method achieves its performance gains without any loss of accuracy to within the local error controlled by the ODE solver. Such speedup allows higher fidelity mechanism chemistry to be coupled with multi-dimensional fluid dynamics simulations.

The density matrix renormalization group algorithm on kilo-processor architectures: Implementation and trade-offs

In the numerical analysis of strongly correlated quantum lattice models one of the leading algorithms developed to balance the size of the effective Hilbert space and the accuracy of the simulation is the density matrix renormalization group (DMRG) algorithm, in which the run-time is dominated by the iterative diagonalization of the Hamilton operator. As the most time-dominant step of the diagonalization can be expressed as a list of dense matrix operations, the DMRG is an appealing candidate to fully utilize the computing power residing in novel kilo-processor architectures.In the paper a smart hybrid CPU–GPU implementation is presented, which exploits the power of both CPU and GPU and tolerates problems exceeding the GPU memory size. Furthermore, a new CUDA kernel has been designed for asymmetric matrix–vector multiplication to accelerate the rest of the diagonalization. Besides the evaluation of the GPU implementation, the practical limits of an FPGA implementation are also discussed.

Solving an elliptic PDE eigenvalue problem via automated multi-level substructuring and hierarchical matrices

We propose a new method for the solution of discretised elliptic PDE eigenvalue problems. The new method combines ideas of domain decomposition, as in the automated multi-level substructuring (short AMLS), with the concept of hierarchical matrices (short $$\fancyscript{H}$$H-matrices) in order to obtain a solver that scales almost optimal in the size of the discrete space. Whereas the AMLS method is very effective for PDEs posed in two dimensions, it is getting very expensive in the three-dimensional case, due to the fact that the interface coupling in the domain decomposition requires dense matrix operations. We resolve this problem by use of data-sparse hierarchical matrices. In addition to the discretisation error our new approach involves a projection error due to AMLS and an arithmetic error due to $$\fancyscript{H}$$H-matrix approximation. A suitable choice of parameters to balance these errors is investigated in examples.

The performance of GRAPE-DR for dense matrix operations

Abstract We describe the implementation and performance of dense matrix multiplication and LU decomposition on the GRAPE-DR SIMD accelerator board. A GRAPE-DR card, with 4 GRAPE-DR chips, has the theoretical peak DP performance of 819 Gflops. Each GRAPE-DR chip has 512 processing elements and operates with 400 MHz clock cycle. each PE can perform one addition and one multiplication in every two clock cycles. The measured performance of matrix multiplication is 730 Gflops for the multiplication of matrices with size 51200 by 2048 and 2048 by 51200. The performance of LU decomposition is 480 Gflops for the problem size of 51200.

A High Performance Multifrontal Code for Linear Solution of Structures Using Multi-Core Microprocessors

A multifrontal code is introduced for the efficient solution of the linear system of equations arising from the analysis of structures. The factorization phase is reduced into a series of interleaved element assembly and dense matrix operations for which the BLAS3 kernels are used. A similar approach is generalized for the forward and back substitution phases for the efficient solution of structures having multiple load conditions. The program performs all assembly and solution steps in parallel. Examples are presented which demonstrate the code's performance on single and dual core processor computers.

An MIMD strategy for quantum mechanical reactive scattering calculations

This article describes a multiprocessor parallel computer program for quantum mechanical reactive scattering calculations by the linear algebraic outgoing wave variational principle. Two aspects of the parallelization are emphasized: (i) domain-based parallelization of multidimensional integrals and (ii) multithreaded and message passing strategies for parallelizing large, dense matrix operations.