Search NASA⌕ Search

SEARCH · Search NASA

Results for “Sparse Matrix Multiplication”

Search indexed NASA NTRS and DOE OSTI research on propulsion, heat transfer, battery materials and energy systems. Follow report and document links to the original sources.

Quote a phrase for an exact phrase match. Source license links do not imply unrestricted reuse.

108 records · Page 6

A supernodal all-pairs shortest path algorithm

We show how to exploit graph sparsity in the Floyd-Warshall algorithm for the all-pairs shortest path (Apsp) problem. Floyd-Warshall is an attractive choice for Apsp on high-performing systems due to its structural similarity to solving dense linear systems and matrix multiplication. However, if sparsity of the input graph is not properly exploited, Floyd-Warshall will perform unnecessary asymptotic work and thus may not be a suitable choice for many input graphs. To overcome this limitation, the key idea in our approach is to use the known algebraic relationship between Floyd-Warshall and Gaussian elimination, and import several algorithmic techniques from sparse Cholesky factorization, namely, fill-in reducing ordering, symbolic analysis, supernodal traversal, and elimination tree parallelism. When combined, these techniques reduce computation, improve locality and enhance parallelism. We implement these ideas in an efficient shared memory parallel prototype that is orders of magnitude faster than an efficient multi-threaded baseline Floyd-Warshall that does not exploit sparsity. Our experiments suggest that the Floyd-Warshall algorithm can compete with Dijkstra's algorithm (the algorithmic core of Johnson's algorithm) for several classes sparse graphs.

Sao, Piyush↗

Extending the QUDA Library with the eigCG Solver

While the incremental eigCG algorithm [ 1 ] is included in many LQCD software packages, its realization on GPU micro-architectures was still missing. In this session we report our experi- ence of the eigCG implementation in the QUDA library. In particular, we will focus on how to employ the mixed precision technique to accelerate solutions of large sparse linear systems with multiple right-hand sides on GPUs. Although application of mixed precision techniques is a well-known optimization approach for linear solvers, its utilization for the eigenvector com- puting within eigCG requires special consideration. We will discuss implementation aspects of the mixed precision deflation and illustrate its numerical behavior on the example of the Wilson twisted mass fermion matrix inversions

Strelchenko, Alexei↗

MATAR: A performance portability and productivity implementation of data-oriented design with Kokkos

There is a need for simple, fast, and memory-efficient multidimensional data structures for dense and sparse storage that arise with numerical methods and in software applications. The data structures must perform equally well across multiple computer architectures, including CPUs and GPUs. For this purpose, we developed MATAR, a C++ software library that allows for simple creation and use of intricate data structures that is also portable across disparate architectures using Kokkos. Here, the performance aspect is achieved by forcing contiguous memory layout (or as close to contiguous as possible) for multidimensional and multi-size dense or sparse MATrix and ARray (hence, MATAR) types. Our results show that MATAR has the capability to improve memory utilization, performance, and programmer productivity in scientific computing. This is achieved by fitting more work into the available memory, minimizing memory loads required, and by loading memory in the most efficient order. This document describes the purpose of the work, the implementation of each of the data types, and the resulting performance both in some simple baseline test cases and in an application code.

97 MATHEMATICS AND COMPUTING↗

Multi-GPU implementation of a VMAT treatment plan optimization algorithm

Purpose: Volumetric modulated arc therapy (VMAT) optimization is a computationally challenging problem due to its large data size, high degrees of freedom, and many hardware constraints. High-performance graphics processing units (GPUs) have been used to speed up the computations. However, GPU’s relatively small memory size cannot handle cases with a large dose-deposition coefficient (DDC) matrix in cases of, e.g., those with a large target size, multiple targets, multiple arcs, and/or small beamlet size. The main purpose of this paper is to report an implementation of a column-generation-based VMAT algorithm, previously developed in the authors’ group, on a multi-GPU platform to solve the memory limitation problem. While the column-generation-based VMAT algorithm has been previously developed, the GPU implementation details have not been reported. Hence, another purpose is to present detailed techniques employed for GPU implementation. The authors also would like to utilize this particular problem as an example problem to study the feasibility of using a multi-GPU platform to solve large-scale problems in medical physics. Methods: The column-generation approach generates VMAT apertures sequentially by solving a pricing problem (PP) and a master problem (MP) iteratively. In the authors’ method, the sparse DDC matrix is first stored on a CPU in coordinate list format (COO). On the GPU side, this matrix is split into four submatrices according to beam angles, which are stored on four GPUs in compressed sparse row format. Computation of beamlet price, the first step in PP, is accomplished using multi-GPUs. A fast inter-GPU data transfer scheme is accomplished using peer-to-peer access. The remaining steps of PP and MP problems are implemented on CPU or a single GPU due to their modest problem scale and computational loads. Barzilai and Borwein algorithm with a subspace step scheme is adopted here to solve the MP problem. A head and neck (H and N) cancer case is then used to validate the authors’ method. The authors also compare their multi-GPU implementation with three different single GPU implementation strategies, i.e., truncating DDC matrix (S1), repeatedly transferring DDC matrix between CPU and GPU (S2), and porting computations involving DDC matrix to CPU (S3), in terms of both plan quality and computational efficiency. Two more H and N patient cases and three prostate cases are used to demonstrate the advantages of the authors’ method. Results: The authors’ multi-GPU implementation can finish the optimization process within ∼1 min for the H and N patient case. S1 leads to an inferior plan quality although its total time was 10 s shorter than the multi-GPU implementation due to the reduced matrix size. S2 and S3 yield the same plan quality as the multi-GPU implementation but take ∼4 and ∼6 min, respectively. High computational efficiency was consistently achieved for the other five patient cases tested, with VMAT plans of clinically acceptable quality obtained within 23–46 s. Conversely, to obtain clinically comparable or acceptable plans for all six of these VMAT cases that the authors have tested in this paper, the optimization time needed in a commercial TPS system on CPU was found to be in an order of several minutes. Conclusions: The results demonstrate that the multi-GPU implementation of the authors’ column-generation-based VMAT optimization can handle the large-scale VMAT optimization problem efficiently without sacrificing plan quality. The authors’ study may serve as an example to shed some light on other large-scale medical physics problems that require multi-GPU techniques.

60 APPLIED LIFE SCIENCES↗

MAGMA: Enabling exascale performance with accelerated BLAS and LAPACK for diverse GPU architectures

MAGMA (Matrix Algebra for GPU and Multicore Architectures) is a pivotal open-source library in the landscape of GPU-enabled dense and sparse linear algebra computations. With a repertoire of approximately 750 numerical routines across four precisions, MAGMA is deeply ingrained in the DOE software stack, playing a crucial role in high-performance computing. Notable projects such as ExaConstit, HiOP, MARBL, and STRUMPACK, among others, directly harness the capabilities of MAGMA. In addition, the MAGMA development team has been acknowledged multiple times for contributing to the vendors’ numerical software stacks. Looking back over the time of the Exascale Computing Project (ECP), we highlight how MAGMA has adapted to recent changes in modern HPC systems, especially the growing gap between CPU and GPU compute capabilities, as well as the introduction of low precision arithmetic in modern GPUs. We also describe MAGMA’s direct impact on several ECP projects. Maintaining portable performance across NVIDIA and AMD GPUs, and with current efforts toward supporting Intel GPUs, MAGMA ensures its adaptability and relevance in the ever-evolving landscape of GPU architectures.

97 MATHEMATICS AND COMPUTING↗

Structural response reconstruction using a system-equivalent singular vector basis

Here, this paper develops a novel method for reconstructing the full-field response of structural dynamic systems using sparse measurements. The singular value decomposition is applied to a frequency response matrix relating the structural response to physical loads, base motion, or modal loads. The left singular vectors form a non-physical reduced basis that can be used for response reconstruction with far fewer sensors than existing methods. The contributions of the singular vectors to measured response are termed singular-vector loads (SVLs) and are used in a regularized Bayesian framework to generate full-field response estimates and confidence intervals. The reconstruction framework is applicable to the estimation of single data records and power spectral densities from multiple records. Reconstruction is successfully performed in configurations where the number of SVLs to identify is less than, equal to, and greater than the number of sensors used for reconstruction. In a simulation featuring a seismically excited shear structure, SVL reconstruction significantly outperforms modal FRF-based reconstruction and successfully estimates full-field responses with as few as two uniaxial accelerometers. SVL reconstruction is further verified in a simulation featuring an acoustically excited cylinder. Finally, response reconstruction and uncertainty quantification are performed on an experimental structure with three shaker inputs and 27 triaxial accelerometer outputs.

42 ENGINEERING↗

Population genomic response to geographic gradients by widespread and endemic fishes of the Arabian Peninsula

Genetic structure within marine species may be driven by local adaptation to their environment, or alternatively by historical processes, such as geographic isolation. The gulfs and seas bordering the Arabian Peninsula offer an ideal setting to examine connectivity patterns in coral reef fishes with respect to environmental gradients and vicariance. The Red Sea is characterized by a unique marine fauna, historical periods of desiccation and isolation, as well as environmental gradients in salinity, temperature, and primary productivity that vary both by latitude and by season. The adjacent Arabian Sea is characterized by a sharper environmental gradient, ranging from extensive coral cover and warm temperatures in the southwest, to sparse coral cover, cooler temperatures, and seasonal upwelling in the northeast. Reef fish, however, are not confined to these seas, with some Red Sea fishes extending varying distances into the northern Arabian Sea, while their pelagic larvae are presumably capable of much greater dispersal. These species must therefore cope with a diversity of conditions that invoke the possibility of steep clines in natural selection. Here, we test for genetic structure in two widespread reef fish species (a butterflyfish and surgeonfish) and eight range-restricted butterflyfishes across the Red Sea and Arabian Sea using genome-wide single nucleotide polymorphisms. We performed multiple matrix regression with randomization analyses on genetic distances for all species, as well as reconstructed scenarios for population subdivision in the species with signatures of isolation. We found that (a) widespread species displayed more genetic subdivision than regional endemics and (b) this genetic structure was not correlated with contemporary environmental parameters but instead may reflect historical events. We propose that the endemic species may be adapted to a diversity of local conditions, but the widespread species are instead subject to ecological filtering where different combinations of genotypes persist under divergent ecological regimes.

59 BASIC BIOLOGICAL SCIENCES↗

GSoFa: Scalable Sparse Symbolic LU Factorization on GPUs

Decomposing a matrix $\mathbf {A}$ into a lower matrix $\mathbf {L}$ and an upper matrix $\mathbf {U}$, which is also known as LU decomposition, is an essential operation in numerical linear algebra. For a sparse matrix, LU decomposition often introduces more nonzero entries in the $\mathbf {L}$ and $\mathbf {U}$ factors than in the original matrix. A symbolic factorization step is needed to identify the nonzero structures of $\mathbf {L}$ and $\mathbf {U}$ matrices. Attracted by the enormous potentials of the Graphics Processing Units (GPUs), an array of efforts have surged to deploy various LU factorization steps except for the symbolic factorization, to the best of our knowledge, on GPUs. This article introduces gSoFa, the first GPU-based symbolic factorization design with the following three optimizations to enable scalable LU symbolic factorization for nonsymmetric pattern sparse matrices on GPUs. First, here we introduce a novel fine-grained parallel symbolic factorization algorithm that is well suited for the Single Instruction Multiple Thread (SIMT) architecture of GPUs. Second, we tailor supernode detection into a SIMT friendly process and strive to balance the workload, minimize the communication and saturate the GPU computing resources during supernode detection. Third, we introduce a three-pronged optimization to reduce the excessive space consumption problem faced by multi-source concurrent symbolic factorization. Taken together, gSoFa achieves up to 31× speedup from 1 to 44 Summit nodes (6 to 264 GPUs) and outperforms the state-of-the-art CPU project, on average, by 5×. Notably, gSoFa also achieves up to 47 percent of the peak memory throughput of a V100 GPU in the Summit Supercomputer.

97 MATHEMATICS AND COMPUTING↗

Educating local radial basis functions using the highest gradient of interest in three dimensional geometries: EDUCATING LOCALIZED RBF-SCHEMES

We present a novel methodology to effectively localize radial basis function approximation methods in three dimensions. The local scheme requires shape parameter-dependent functions that can be used to approximate gradients of scattered data and to solve partial differential operators. The optimum shape parameter is obtained from the highest gradient of interest, where a known analytical function, when boundary conditions are not present, or a shape parameter-free global approximation are used to educate the localized scheme. The later option is applicable to problems where the operator needs to be solved multiple times, like in time evolution or stochastic integration. Past shape parameter's optimizations, for two-dimensional domains, based on the condition number of the interpolant matrix, were unable to provide satisfactory approximations. The applicability of our method is illustrated in the context of an analytical expression interpolation and during a Ginzburg-Landau relaxation of a free energy functional. In general, the optimum shape parameter depends on geometry, node distribution, and density, whereas the approximation errors decrease as the node density and the local stencil size increase. The effective localization of radial basis functions motivates its use in moving boundary problems and accelerates solutions through sparse matrix solvers.

free energy functional↗

Comparison of sparse domain approaches for 4D SPECT dynamic image reconstruction

Purpose: Dynamic imaging (DI) provides additional diagnostic information in emission tomography in comparison to conventional static imaging at the cost of being computationally more challenging. Dynamic single photon emission computed tomography (SPECT) reconstruction is particularly difficult because of the limitations in the sampling geometry present in most existing scanners. We have developed an algorithm Spline Initialized Factor Analysis of Dynamic Structures (SIFADS) that is a matrix factorization method for reconstructing the dynamics of tracers in tissues and blood directly from the projections in dynamic cardiac SPECT, without first resorting to any 3D reconstruction. Methods: SIFADS is different from "pure" factor analysis in dynamic structures (FADS) in that it employs a dedicated spline-based pre-initialization. In this paper, we analyze the convergence properties of SIFADS and FADS using multiple metrics. The performances of the two approaches are evaluated for numerically simulated data and for real dynamic SPECT data from canine and human subjects. Results: For SIFADS, metrics analyzed for reconstruction algorithm convergence show better features of the metric curves vs iterations. In addition, SIAFDS provides better tissue segmentations than that from pure FADS. Measured computational times are also typically better for SIFADS implementations than those with pure FADS. Conclusion: The analysis supports the utility of the pre-initialization of a factorization algorithm for better dynamic SPECT image reconstruction.

46 INSTRUMENTATION RELATED TO NUCLEAR SCIENCE AND ↗

HPC-NMF: A High-Performance Parallel Algorithm for Nonnegative Matrix Factorization

NMF is a useful tool for many applications in different domains such as topic modeling in text mining, background separation in video analysis, and community detection in social networks. Despite its popularity in the data mining community, there is a lack of efficient distributed algorithms to solve the problem for big data sets. We propose a high-performance distributed-memory parallel algorithm that computes the factorization by iteratively solving alternating non-negative least squares (NLS) subproblems for $\WW$ and $\HH$. It maintains the data and factor matrices in memory (distributed across processors), uses MPI for interprocessor communication, and, in the dense case, provably minimizes communication costs (under mild assumptions). As opposed to previous implementation, our algorithm is also flexible: It performs well for both dense and sparse matrices, and allows the user to choose any one of the multiple algorithms for solving the updates to low rank factors $\WW$ and $\HH$ within the alternating iterations.

Kannan, Ramakrishnan↗

Parallel implementation and performance optimization of the configuration-interaction method

The configuration-interaction (CI) method, long a popular approach to describe quantum many-body systems, is cast as a very large sparse matrix eigenpair problem with matrices whose dimension can exceed one billion. Such formulations place high demands on memory capacity and memory bandwidth - - two quantities at a premium today. In this paper, we describe an efficient, scalable implementation, BIGSTICK, which, by factorizing both the basis and the interaction into two levels, can reconstruct the nonzero matrix elements on the fly, reduce the memory requirements by one or two orders of magnitude, and enable researchers to trade reduced resources for increased computational time. We optimize BIGSTICK on two leading HPC platforms - - the Cray XC30 and the IBM Blue Gene/Q. Specifically, we not only develop an empirically-driven load balancing strategy that can evenly distribute the matrix-vector multiplication across 256K threads, we also developed techniques that improve the performance of the Lanczos reorthogonalization. Combined, these optimizations improved performance by 1.3-8× depending on platform and configuration.

97 MATHEMATICS AND COMPUTING↗

A high-performance parallel algorithm for nonnegative matrix factorization

Non-negative matrix factorization (NMF) is the problem of determining two non-negative low rank factors W and H, for the given input matrix A, such that A ≈ WH. NMF is a useful tool for many applications in different domains such as topic modeling in text mining, background separation in video analysis, and community detection in social networks. Despite its popularity in the data mining community, there is a lack of efficient distributed algorithms to solve the problem for big data sets. We propose a high-performance distributed-memory parallel algorithm that computes the factorization by iteratively solving alternating non-negative least squares (NLS) subproblems for W and H. It maintains the data and factor matrices in memory (distributed across processors), uses MPI for interprocessor communication, and, in the dense case, provably minimizes communication costs (under mild assumptions). As opposed to previous implementations, our algorithm is also flexible: (1) it performs well for both dense and sparse matrices, and (2) it allows the user to choose any one of the multiple algorithms for solving the updates to low rank factors W and H within the alternating iterations. We demonstrate the scalability of our algorithm and compare it with baseline implementations, showing significant performance improvements.

Kannan, Ramakrishnan↗

Mathematical foundations of the GraphBLAS

The GraphBLAS standard (GraphBlas.org) is being developed to bring the potential of matrix-based graph algorithms to the broadest possible audience. Mathematically, the GraphBLAS defines a core set of matrix-based graph operations that can be used to implement a wide class of graph algorithms in a wide range of programming environments. This study provides an introduction to the mathematics of the GraphBLAS. Graphs represent connections between vertices with edges. Matrices can represent a wide range of graphs using adjacency matrices or incidence matrices. Adjacency matrices are often easier to analyze while incidence matrices are often better for representing data. Fortunately, the two are easily connected by matrix multiplication. A key feature of matrix mathematics is that a very small number of matrix operations can be used to manipulate a very wide range of graphs. This composability of a small number of operations is the foundation of the GraphBLAS. A standard such as the GraphBLAS can only be effective if it has low performance overhead. Finally, performance measurements of prototype GraphBLAS implementations indicate that the overhead is low.

97 MATHEMATICS AND COMPUTING↗

Optimizing shift selection in multilevel Monte Carlo for disconnected diagrams in lattice QCD

The calculation of disconnected diagram contributions to physical signals is a computationally expensive task in Lattice QCD. To extract the physical signal, the trace of the inverse Lattice Dirac operator, a large sparse matrix, must be stochastically estimated. Because the variance of the stochastic estimator is typically large, variance reduction techniques must be employed. Multilevel Monte Carlo (MLMC) methods reduce the variance of the trace estimator by utilizing a telescoping sequence of estimators. Frequency Splitting is one such method that uses a sequence of inverses of shifted operators to estimate the trace of the inverse lattice Dirac operator, however there is no a priori way to select the shifts that minimize the cost of the multilevel trace estimation. Here we present a sampling and interpolation scheme that is able to predict the variances associated with Frequency Splitting under displacements of the underlying space time lattice. The interpolation scheme is able to predict the variances to high accuracy and therefore chooses shifts that correspond to an approximate minimum of the cost for the trace estimation. We show that Frequency Splitting with the chosen shifts displays significant speedups over multigrid deflation, and that these shifts can be used for multiple configurations within the same ensemble with no penalty to performance.

97 MATHEMATICS AND COMPUTING↗

MPI-FAUN: An MPI-Based Framework for Alternating-Updating Nonnegative Matrix Factorization

Non-negative matrix factorization (NMF) is the problem of determining two non-negative low rank factors W and H, for the given input matrix A, such that A≈WH. NMF is a useful tool for many applications in different domains such as topic modeling in text mining, background separation in video analysis, and community detection in social networks. Despite its popularity in the data mining community, there is a lack of efficient parallel algorithms to solve the problem for big data sets. The main contribution of this work is a new, high-performance parallel computational framework for a broad class of NMF algorithms that iteratively solves alternating non-negative least squares (NLS) subproblems for W and H. It maintains the data and factor matrices in memory (distributed across processors), uses MPI for interprocessor communication, and, in the dense case, provably minimizes communication costs (under mild assumptions). The framework is flexible and able to leverage a variety of NMF and NLS algorithms, including Multiplicative Update, Hierarchical Alternating Least Squares, and Block Principal Pivoting. Our implementation allows us to benchmark and compare different algorithms on massive dense and sparse data matrices of size that spans from few hundreds of millions to billions. We demonstrate the scalability of our algorithm and compare it with baseline implementations, showing significant performance improvements. The code and the datasets used for conducting the experiments are available online.

2D↗

Analog Systems for Edge Optimization

Over the past decade, analog computing has the subject of substantial research interest providing a path toward improved computational efficiency in the post-Dennard era. Analog matrix vector multiplication (MVM) accelerators provide a popular approach given the ubiquity of MVM operations in numerous applications. However, historically analog computing systems can struggle with applications requiring high precision due to the inherent susceptibility of these systems to analog non-idealities. Therefore, prior work on analog systems has focused either on applications known to be tolerant of limited precision (e.g., neural network inference), or using expensive techniques to emulate high-precision using many analog MVM operations. In this work, we propose an alternative approach. Motivated by recent advances in inexact nonlinear solvers and optimizers, we explore the potential of co-designing optimization algorithms which can take full advantage of the fundamentally inexact analog MVM operations. To enable these co-designed algorithms we also develop a general mathematical theory of the precision and energy efficiency of analog operations, and a new system architecture for tightly-coupled analog and digital computation. Finally, we examine the applicability of analog computing to a wider class of symmetric positive definite systems and find potential in using analog operations as a sparse approximate inverse preconditioner. With these core innovations, this project provides a path toward effectively implementing optimization algorithms on power-constrained autonomous and semi-autonomous systems.

97 MATHEMATICS AND COMPUTING↗

Sparse Control Synthesis for Uncertain Responsive Loads With Stochastic Stability Guarantees

In this report, recent studies have demonstrated the potential of flexible loads in providing frequency response services, predominantly due to their availability and cost-effectiveness. However, uncertainty and variability in various weather-related and end-use behavioral factors often impact the reliability of demand-side control performance. This work addresses this problem with the design of a demand-side control to achieve frequency response under load uncertainties. Our approach involves modeling the load uncertainties via stochastic processes that appear as both multiplicative and additive in the power system dynamics. Recently developed mean square exponential stability (MSES) results for continuous-time linear stochastic systems are applied to pose the control synthesis problem which results in an LMI-based optimization problem. Additional costs and constraints are added to the LMI-based controller synthesis to ensure MSES, improve closed-loop transient performance, maximize tolerable uncertainties, and promote sparsity in the controller. Additionally, the fundamental limitations between the tolerable uncertainties and control efforts while ensuring MSES are discussed. Further, the control synthesis problem for the case of the full-state measurement is generalized to the case of partial-state measurements. The proposed control synthesis is illustrated on an IEEE 39 bus system with rigorous studies to demonstrate the role of sparsity, closed-loop transient performance, tolerable uncertainties, and control efforts while ensuring MSES and achieving frequency response.

42 ENGINEERING↗