SMP Based Solver for Large Binary Systems

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.

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.

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.

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.

A parallel electrostatic solver for the VORPAL code

A projection method for solving nonsymmetric linear systems on multiprocessors

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.

Algorithms for the sparse matrix-vector multiplication (shortly SpM×V ) are important building blocks in solvers of sparse systems of linear equations. Due to matrix sparsity, the memory access patterns are irregular and the utilization of a cache suffers from low spatial and temporal locality. To reduce this effect, the register blocking formats were designed. This paper introduces a new combined format, for storing sparse matrices that extends possibilities of the variable-sized register blocking format.

The polarizable continuum model (PCM) is a well established method for computing solvation effects. Its attractiveness comes from the fact that the solution is modelled as a continuum instead of modelling each atom individually. This allows for more complex simulations. Nevertheless, there are still a couple of challenges to face. One of them being the computational time required at the limit of large systems, or when high accuracy is needed. The wavelet boundary element method can be used to overcome these problems, provided a reliable cavity generator is used. The PCM calculates the molecular free energy in solution as the sum of electrostatic, dispersion-repulsion and the cavitation energy by solving the underlying partial differential equations. Those differential equations can be transformed into integral equations, solely defined on the boundary of the molecular cavity, by applying the integral equation formalism (IEFPCM). The integral equations can then be discretized using the wavelet boundary element method. The resulting sparse system of linear equations can be solved reliably by iterative solvers, leading to high accuracy solutions even for large systems. The challenge of this approach lies in the generation and refinement of the molecular cavities on which the interactions take place. Since wavelets are defined as a tensor product of one-dimensional functions, the wavelet boundary element method needs a discretization into quadrangular patches of the molecular cavity. This feature is missing in common commercial mesh generation tools, which focus on generating triangular meshes. Making use of characteristic functions for defining which points are inside the cavity, a model for generating the molecular surface independently of its exact description is achieved. We present here a way for generating triangular meshes which are subsequently merged into quadrangular ones. We apply geometrical quantities that measure the fitness of the shapes involved and improve the position of individual points iteratively. Thereby, we generate quadrangular meshes with an unprecedented quality, which ultimately can be refined in a hierarchical manner for the use of wavelets. Overall, this method can generate reliable parametrizations on a variety of smooth cavities. With the help of the wavelet boundary element method, the PCM equations can thus be solved with high accuracy also on large systems.

To analyze the behavior of semiconductor devices numerically, we need to solve extremely large and sparse linear systems of equations as fast as possible. In this paper, parallel processing architectures using systolic arrays are proposed as a rapid solver of linear equations. The architectures are designed to reduce the number of processing elements, utilizing the regularity and the sparseness of coefficient matrices.The proposed systolic arrays consist of sixty‐four cells for incomplete matrix decomposition and of seven cells for matrix‐vector multiplication. This number is independent from the number N of grid points discretized for a three‐dimensional device. The processing time of the array to solve a system of equations is proportional to N. Furthermore, the architectures are revised for faster parallel computation. The acceleration is based on the fact that a coefficient matrix is decomposed of smaller matrices. By the acceleration for matrix‐vector multiplication, required cells are 8 × N1/3, but the processing time is proportional to N2/3.

Performance assessment and field-development planning for carbon geo-sequestration relies on reservoir simulation. The standard approach is to use regular (corner-point) grids in space and a fully implicit discretisation in time, solving for all primary unknowns at once using Newton’s method. While the former limits the physical realism of the simulation model built from the geomodel, the latter leads to a large ill-conditioned system of equations that is inefficient to solve. We overcome the regularisation issue via a hierarchy of fully unstructured, geobody-conforming, finite element meshes which can be refined until mesh convergence is achieved. To address the time discretisation issue we employ, for the first time, the linearly implicit extrapolation scheme (LIMEX) to solve the highly non-linear coupled two-phase flow and reactive transport equations. To solve the arising large sparse systems of linear equations, we apply the geometric multigrid (GMG) method that demonstrates optimal, linear complexity and allows an efficient parallelization on supercomputers. Another novel feature of our formulation is the consideration of the kinetics of mass transfer between the carbonic and aqueous phases. This approach removes the first-order dependence on mesh refinement of CO2 that is dissolved, a characteristic feature of standard equilibrium models. We demonstrate the parallel scalability of our simulation framework with an implementation based on the UG4 platform. Proof-of-concept results accurately capture key features of CO2 migration including filtration by capillary barriers and convective dissolution of CO2 at the base of the plume.

We consider matrices arising out of the one-way dissection method for solving large sparse systems of linear equations. The systems that we consider are those that may have singular diagonal blocks. Such systems arise in certain fluid flow problems. Gunzberger and Nicholaides proposed a method for resolving the singularity in the diagonal blocks. This method uses the Moore-Penrose pseudoinverse. We propose two improvements to the Gunzberger-Nicholaides procedure: the substitution of a weighted pseudoinverse for the Moore-Penrose pseudoinverse; a more elegant implementation of the back substitution procedure. A stability analysis of both our procedure and the Gunzberger-Nicholaides procedure is given. Both our analysis and empirical tests show that our method has better numerical stability properties than the Gunzberger-Nicholaides procedure. We also implement our procedure on Intel iPSC/1 Hypercube. Our improvement to the back substitution method makes the natural parallelism in the problem easier to exploit. 13 refs., 2 tabs.

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.