Sparse Matrix-Vector Multiplication - Final Solution?

During the past few years, algorithmic improvements alone have reduced the time required for the direct solution of unsymmetric sparse systems of linear equations by almost an order of magnitude. This paper compares the performance of some well-known software packages for solving general sparse systems. In particular, it demonstrates the consistently high level of performance achieved by WSMP---the most recent of such solvers. It compares the various algorithmic components of these solvers and discusses their impact on solver performance. Our experiments show that the algorithmic choices made in WSMP enable it to run more than twice as fast as the best among similar solvers and that WSMP can factor some of the largest sparse matrices available from real applications in only a few seconds on a 4-CPU workstation. Thus, the combination of advances in hardware and algorithms makes it possible to solve those general sparse linear systems quickly and easily that might have been considered too large until recently.

This paper presents a novel hardware solution for solving sparse matrices using emulation technology. The proposed solution introduces a bio-inspired technique, which is Genetic Algorithm (GA) to solve the sparse system of linear equations (SLEs). The proposed solution can be used with any physical system that can be modeled as SLEs. Computations involved in Finite Element Method (FEM) consume too much time, which affects the final time-to-market value. Profiling shows that the most time-consuming part in the simulation process is the solver part which is responsible for solving the resultant SLEs generated from the FEM. The total number of equations may reach thousands or millions of linear equations. Hardware implementation of GA and making a good use of parallelism in implementing the architecture accelerate the time taken by the solver part. Solving the system by classical numerical methods, such as conjugate gradient (CG) implemented in software is performed sequentially, here the architecture allows the parallelism in performing GA operators, such as selection, crossover, and mutation. The parallelism in the proposed architecture is the superiority point which we rely on. The results show that our proposed method to solve the SLEs outperforms the classical CG software implementation on run time and number of equations.

Solving a large sparse system of linear equations (SSLE)over the binary field is an important step in the general number field seive (GNFS) algorithm that is often used to factorize large integers used in the RSA algorithm. Block Lanczos method proposed by Montgomery is an efficient and reliable method for solving such systems. A number of implementations of Block Lanczos method on distributed systems have been proposed in past. This paper discusses the parallel Montgomery’s Block Lanczos method for solving binary system of equations on shared multiprocessors(SMP). A simple experiment shows that the speed of convergence of this algorithm is dependent on a wise choice of initial guess for kick-starting the algorithm. A somewhat dense choice of the initial guess is suggested rather than a random choice for faster convergence. Devoid of any communication overheads of a loosely coupled system, the improved method gives good performance on SMP systems which can provide an able alternative to the more popular distributed systems. Details of implementation of this algorithm on SMP and experimental results are also provided in the paper.

Computer subroutines for solving large sparse systems of linear equations usually involve fairly complicated data structures and storage management. In most cases, the user of such subroutines simply wants top solve his problem and does not have to understand how the storage management is done, or how the matrix components are stored. Thus software for solving sparse systems of linear equations should be designed in such a way that it can insulate the user from these considerations, while still allowing the software to be used in a variety of ways.

During my visit to LLNL during the summer of 2010, I worked on algebraic multilevel solvers for large sparse systems of linear equations arising from discretizations of partial differential equations. The particular solver of interest is based on ILU decomposition. The setup phase for this AMG solve is just the single ILU decomposition, and its corresponding error matrix. Because the ILU uses a minimum degree or similar sparse matrix ordering, most of the fill-in, and hence most of the error, is concentrated in the lower right corner of the factored matrix. All of the major multigrid components - the smoother, the coarse level correction matrices, and the fine-to-coarse and coarse-to-fine rectangular transfer matrices, are defined in terms of various blocks of the ILU factorization. Although such a strategy is not likely to be optimal in terms of convergence properties, it has a relatively low setup cost, and therefore is useful in situations where setup costs for more traditional AMG approaches cannot be amortized over the solution of many linear systems using the same matrix. Such a situation arises in adaptive methods, where often just one linear system is solved at each step of an adaptive feedback loop, or in solving nonlinear equations by approximate Newton methods, where the approximate Jacobian might change substantially from iteration to iteration. In general terms, coarse levels are defined in terms of successively smaller lower right blocks of the matrix, typically decreasing geometrically in order. The most difficult issue was the coarse grid correction matrix. The preconditioner/smoother for a given level is just the corresponding lower right blocks of the ILU factorization. The coarse level matrix itself is just the Schur complement; this matrix is not known exactly using just the ILU decomposition in the setup phase. Thus we approximate this matrix using various combinations of the preconditioning matrix and the error matrix. During my visit, several approximations of this type were implemented and tested. While some improved the convergence rate of the overall method, these gains had to be balanced against the additional costs involved in creating and applying these matrices. By this more stringent criterion, none of the improved approximations could be characterized as an unqualified success.

The DEFLATED-GMRES(m, k) method is one of the major iterative solvers for the large sparse linear systems of equations, x = b. This algorithm assembles a preconditioner adaptively for the GMRES(m) method based on eigencomponents gathered from the Arnoldi process during iterations. It is usually known that if a restarted GMRES(m) method is used to solve linear systems of equations, the information of the smallest eigencomponents is lost at each restart and the super-linear convergence may also be lost. In this paper, we propose an adaptive procedure that combines the DEFLATED-GMRES(m, k) algorithm and the determination of a restart frequency automatically. It is shown that a new algorithm combining elements of both will reduce the negative effects of the restarted procedure. The numerical experiments are presented on three test problems by using the MIMD parallel machine AP3000. From these numerical results, we show that the proposed algorithm leads to faster convergence than the conventional DEFLATED-GMRES(m, k) method. Copyright © 2000 John Wiley & Sons, Ltd.

Data locality exists everywhere in the memory hierarchy. Most applications show temporal and spatial locality, and computer system and architecture designers utilize this property to improve the system performance with better data layout, prefetching, and scheduling. The locality property can be represented by memory access patterns, which records the time and frequency of accessed addresses. From the security perspective, if an attacker can trace the access pattern, sensitive information inside of the application could be observed and leaked. Oblivious RAM is one of the most effective solutions to mitigate the access pattern leakage on the system, which adds redundant data blocks in space and time. With ORAM protection, the intrinsic data locality is broken by the randomly stored data. Therefore, the application cannot gain any performance benefits from locality if the ORAM protocol is used. In this work, we would like to study the potential to support multi-range accesses with new storage and access efficient ORAM construction. Our proposed designs include two major schemes: Lite-rORAM, which minimize the storage overhead of existing rORAM; and Hybrid-rORAM, which support multiple ranges accesses with minimum storage overhead. We achieve the goal to preserve the locality for consecutive data blocks with different ranges in the application while obfuscates the access pattern as well. We tested our proposed schemes with different workloads on local and remote backends. The experimental results show that, in the best case, our proposed ORAM construction can reduce the data block retrieval time to 0.24x of the baseline Path ORAM, with 87.5% storage overhead reduction compared to rORAM.

It is known that the main computational costs in the FEM implementation are computing the global system of linear equations (material response computation and stiffness matrix ensemble) and its solution. With the aim of speeding up the total solving boundary value problem time with the help of the developed solver, algorithms of the parallel stiffness matrix ensemble and the system of linear equations solution were implemented. In this paper the results of using GPU to accelerate computations related to the solution of direct boundary value problems of crystal plasticity are present. The existing algorithms for solving the sparse systems of linear equations and parallel ensemble of the global stiffness matrix were implemented on both the CPU and GPU. As a result, a total acceleration of up to 10 times was obtained compared to the software implementation on the CPU.

Many applications in science and engineering give rise to large sparse linear systems of equations that need to be solved as efficiently as possible. As the size of the problems of interest increases, it can become necessary to consider exploiting multiprocessors to solve these systems. We report on the design and development of parallel frontal solvers for the numerical solution of large sparse linear systems. Three codes have been developed for the mathematical software library HSL (www.cse.clrc.ac.uk/Activity/HSL). The first is for unsymmetric finite-element problems; the second is for symmetric positive definite finite-element problems; and the third is for highly unsymmetric linear systems such as those that arise in chemical process engineering. In each case, the problem is subdivided into a small number of loosely connected subproblems and a frontal method is then applied to each of the subproblems in parallel. We discuss how our software is designed to achieve the goals of portability, ease of use, efficiency, and flexibility, and illustrate the performance using problems arising from real applications.

Element-by-element iterative solution methods have been shown to be effective for solving large, sparse linear systems of equations with a positive definite stiffness matrix. In this paper the iterative method GMRES will be described. This solution method is also suited for solving nonsymmetric systems of equations. As acceleration technique an element-by-element preconditioner has been employed. Attention is focused mainly on the convergence behaviour of the method. The calculations described in this paper have been performed with the finite element package DIANA.

Evolutionary multimodal optimization using the principle of locality

P aS tiX: a high-performance parallel direct solver for sparse symmetric positive definite systems

Introduction. Kron [1] has developed a technique for tearing large, sparse linear systems of algebraic equations into smaller systems, then putting the solutions of these smaller systems together to form the solution of the original problem. This can result in a great savings over solving the whole system in one piece. The trick is to find where to tear so that on the one hand the remaining pieces are small and easily solvable, but on the other hand not so much violence is done to the equations so that the task of putting their solutions back together becomes too difncult. Kron's pro? cedure for solving the pieces and putting their solutions together is based upon representing systems of equations as electric networks. He makes his choice of where to tear from his insight into these networks and the physical processes the equations represent. In the hands of someone such as Kron who has an understanding of the methods and good insight into the prob? lem, this technique has had outstanding success. The widespread use of Kron's technique has perhaps been limited severely for lack of practitioners skilled in this art of determining effective ways of tearing specific problems. Many problems are too large to submit to such insight. It is the purpose of this paper to contribute to the development of procedures which can be performed on a computer for analyzing the structures of the systems of equations themselves as an aid in choosing how to break them up for easier solution. The methods developed here require only that we consider what variables occur in each equation, but not how they appear. Thus, the tech? niques are not limited to linear systems.

A projection method for solving nonsymmetric linear systems on multiprocessors

Image acquisition of magnetic resonance imaging (MRI) can be accelerated by using multiple receiving coils simultaneously. The problem of reconstructing an unaliased image from partially sampled k-space data can be formulated as a large system of sparse linear equations. The k-space sparse matrix (kSPA) algorithm proposes to solve the system of equations by finding a sparse approximate inverse. This algorithm has been shown to accelerate the image reconstruction for a large number of coils. The original kSPA algorithm requires knowledge of coil sensitivities. Here, we propose and demonstrate an auto-calibrated kSPA algorithm that does not require the explicit computation of the coil sensitivity maps. We have also shown that calibration data, in principle, can be acquired at any region of k-space. This property applies to arbitrary sampling trajectories and all reconstruction algorithms based on k-space. In practice, because of its higher SNR, calibration data acquired at the center of k-space performed more favorably. Such auto-calibration can be advantageous in cases where an accurate sensitivity map is difficult to obtain.