Distributed Matrix-Free Solution of Large Sparse Linear Systems over Finite Fields

New features of our DSC system for distributing a symbolic computation task over a network of processors are described. A new scheduler sends parallel subtasks to those compute nodes that are best suited in handling the added load of CPU usage and memory. Furthermore, a subtask can communicate back to the process that spawned it by a co-routine style calling mechanism. Two large experiments are described in this improved setting. We have implemented an algorithm that can prove a number of more than 1,000 decimal digits prime in about 2 months elapsed time on some 20 computers. A parallel version of a sparse linear system solver is used to compute the solution of sparse linear systems over finite fields. We are able to find the solution of a 100,000 by 100,000 linear system with about 10.3 million non-zero entries over the Galois field with 2 elements using 3 computers in about 54 hours CPU time.

Efficient approximate solution of sparse linear systems

In this dissertation we develop numerical algorithms and software tools to facilitate parallel solutions of nonlinear programming (NLP) problems. In particular, we address large-scale, block-structured problems with an intrinsic decomposable configuration. These problems arise in a great number of engineering applications, including parameter estimation, optimal control, network optimization, and stochastic programming. The structure of these problems can be leveraged by optimization solvers to accelerate solutions and overcome memory limitations, and we propose variants to two classes of optimization algorithms: augmented Lagrangian (AL) schemes and Schur-complement interior-point methods. The convergence properties of augmented Lagrangian decomposition schemes like the alternating direction method of multipliers (ADMM) and progressive hedging (PH) are well established for convex optimization but convergence guarantees in non-convex settings are still poorly understood. In practice, however, ADMM and PH often perform satisfactorily in complex non-convex NLPs. In this work, we study connections between the method of multipliers (MM), ADMM, and PH to derive benchmarking metrics that explain why PH and ADMM work in practice. We illustrate the concepts using challenging dynamic optimization problems. Our exposition seeks to establish more formalism in benchmarking ADMM, PH, and AL schemes and to motivate algorithmic improvements. The effectiveness of nonlinear interior-point solvers for solving large-scale problems relies quite heavily on the solution of the underlying linear algebra systems. The schur-complement decomposition is very effective for parallelizing the solution of linear systems with modest coupling. However, for systems with large number of coupling variables the schur-complement method does not scale favorably. We implement an approach that uses a Krylov solver (GMRES) preconditioned with ADMM to solve block-structured linear systems that arise in the interior-point method. We show that this ADMM-GMRES approach overcomes the well-known scalability issues of Schur decomposition. One important drawback of using decomposition approaches like ADMM and PH is their convergence rate. Unlike Schur-complement interior-point algorithms that have super-linear convergence, augmented Lagrangian approaches typically exhibit linear and sublinear rates. We exploit connections between ADMM and the Schur-complement decomposition to derive an accelerated version of ADMM. Specifically, we study the effectiveness of performing a Newton-Raphson algorithm to compute multiplier estimates for augmented Lagrangian methods. We demonstrate using two-stage stochastic programming problems that our multiplier update achieves convergence in fewer iterations for MM on general nonlinear problems. In the case of ADMM, the newton update significantly reduces the number of subproblem solves for convex quadratic programs (QPs). Moreover, we show that using newton multiplier updates makes the method robust to the selection of the penalty parameter. Traditionally, state-of-the-art optimization solvers are implemented in low-level programming languages. In our experience, the development of decomposition algorithms in these frameworks is challenging. They present a steep learning curve and can slow the development and testing of new numerical algorithms. To mitigate these challenges, we developed PyNumero, a new open source framework implemented in Python and C++. The package seeks to facilitate development of optimization algorithms for large-scale optimization within a high-level programming environment while at the same time minimizing the computational burden of using Python. The efficiency of PyNumero is illustrated by implementing algorithms for problems arising in stochastic programming and optimal control. Timing results are presented for both serial and parallel implementations. Our computational studies demonstrate that with the appropriate balance between compiled code and Python, efficient implementations of optimization algorithms are achievable in these high-level languages.

We study the role of preconditioning strategies recently developed for coercive problems in connection with a two-step iterative method, based on the Hermitian skew-Hermitian splitting (HSS) of the coefficient matrix, proposed by Bai, Golub and Ng for the solution of nonsymmetric linear systems whose real part is coercive. As a model problem we consider Finite Differences (FD) matrix sequences {A n (a,p)} n discretizing the elliptic (convection-diffusion) problem with Ω being a plurirectangle of R d with a(x) being a uniformly positive function and p(x) denoting the Reynolds function: here for plurirectangle we mean a connected union of rectangles in d dimensions with edges parallel to the axes. More precisely, in connection with preconditioned HSS/GMRES like methods, we consider the preconditioning sequence {P n (a)} n , P n (a):= D n 1/2(a)A n (1,0) D n 1/2(a) where D n (a) is the suitably scaled main diagonal of A n (a,0). If a(x) is positive and regular enough, then the preconditioned sequence shows a strong clustering at unity so that the sequence {P n (a)} n turns out to be a superlinear preconditioning sequence for {A n (a,0)} n where A n (a,0) represents a good approximation of Re(A n (a,p)) namely the real part of A n (a,p). The computational interest is due to the fact that the preconditioned HSS method has a convergence behavior depending on the spectral properties of {P n -1(a)Re(A n (a,p))} n ≈ {P n -1(a)A n (a,0)} n : therefore the solution of a linear system with coefficient matrix A n (a,p) is reduced to computations involving diagonals and to the use of fast Poisson solvers for {A n (1,0)} n . Some numerical experimentations confirm the optimality of the discussed proposal and its superiority with respect to existing techniques.

This chapter treats the theory of ordinary differential equations, both linear and nonlinear.Sections 1–4 establish existence and uniqueness theorems for ordinary differential equations. The first section gives some examples of first-order equations, mostly nonlinear, to illustrate certain kinds of behavior of solutions. The second section shows, in the presence of continuity for a vectorvalued F satisfying a “Lipschitz condition,” that the first-order system y′ = F(t, y) has a unique local solution satisfying an initial condition y (t0) = y0. Since higher-order equations can always be reduced to first-order systems, these results address existence and uniqueness for nth-order equations as a special case. Section 3 shows that the solutions to a system depend well on the initial condition and on any parameters that are present in F. Section 4 applies these results to existence of integral curves for a vector field and to construction of coordinate systems from families of integral curves.Sections 5–8 concern linear systems. Section 5 shows that local solutions of linear systems may be extended to global solutions and that in the homogeneous case the vector space of global solutions has dimension equal to the size of the system. The method of variation of parameters reduces the solution of any linear system to the solution of a homogeneous linear system. Sections 6–7 identify explicit solutions to nth-order linear equations and first-order linear systems. The “Jordan canonical form” of a square matrix plays a role in the case of a system. Section 8 discusses power-series solutions to second-order homogeneous linear equations whose coefficients are given by convergent power series, as well as solutions that arise in the case of regular singular points. Two kinds of special functions are mentioned that result from this study-Legendre polynomials and Bessel functions.KeywordsOrdinary Differential EquationUniqueness TheoremDiagonal EntryLipschitz ConditionIntegral CurfThese keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

This paper describes and tests a parallel implementation of a factorized approximate inverse preconditioner (FSAI) to accelerate iterative linear system solvers. Such a preconditioner reveals an efficient accelerator of both Conjugate gradient and BiCGstab iterative methods in the parallel solution of large linear systems arising from the discretization of the advection-diffusion equation. The resulting message passing code allows the solution of large problems leading to a very cost-effective algorithm for the solution of large and sparse linear systems.

A recursive algorithm for generating the transition matrices of multistation multiserver exponential reliable queueing networks

Developments and trends in the parallel solution of linear systems

The goal of this paper is to develop a parallel algorithm for the direct solution of large sparse linear systems and integrate it into domain decomposition methods. The computational effort for these linear systems, often encountered in numerical simulation of structural mechanics problems by finite element codes, is very significant in terms of run-time and memory requirements.In this paper, a two-level parallelism is exploited. The exploitation of the lower level of parallelism is based on the development of a parallel direct solver with a nested dissection algorithm and to introduce it into the FETI methods. This direct solver has the advantage of handling zero-energy modes in floating structures automatically and properly. The upper level of parallelism is a coarse-grain parallelism between substructures of FETI. Some numerical tests are carried out to evaluate the performance of the direct solver.

This work deals with the solution of large non-Hermitian linear systems on desktop workstations with multiple graphics processing units (GPUs). While our implementation is motivated by the need to accelerate volume conductor modeling for bioelectrical brain imaging, the problem itself is common in scientific computing. Whenever a complex partial differential equation is numerically solved, a typically non-Hermitian sparse complex linear system needs to be solved. For problem sizes in the millions, this can take a long time even with highly optimized CPU-based solvers. Our GPU-accelerated solver outperforms an optimized OpenMP-based reference running on two quad-core CPUs by a factor of up to 31× in single precision and up to 7× in double precision, at the cost of a very modest hardware upgrade of two dual-GPU GTX 295 graphics cards. A pair of stronger Fermi GPUs (GTX 480) achieves speedups of 30× in single precision and 15× in double precision.

Accelerated generalized successive overrelaxation method for least squares problems

Black box techniques [12] are enabling exact linear algebra computations of a scale well beyond anything previously possible. The development of new and interesting algorithms has proceeded apace for the past two decades. It is time for the dissemination of these algorithms in an easily used software library so that the mathematical community may readily take advantage of their power. LinBox is that library. In this paper, we describe the design of this generic library, sketch its current range of capabilities, and give several examples of its use. The examples include a solution of Trefethen’s “Hundred Digit Challenge” problem #7 [14] and the computation of all the homology groups of simplicial complexes using the Smith normal form [8]. Exact black box methods are currently successful on sparse matrices with hundreds of thousands of rows and columns and having several million nonzero entries. The main reason large problems can be solved by black box methods is that they require much less memory in general than traditional eliminationbased methods do. This fact is widely used in the numerical computation area. We refer for instance to the templates for linear system solution and eigenvalue problems [2,1]. This has also led the computer algebra community to a considerable interest in black box methods. Since Wiedemann’s seminal paper [16], many developments have been proposed especially to adapt Krylov or Lanczos methods to fast exact algorithms. We refer to [5] and references therein for a review of problems and solutions. LinBox supplies efficient black box solutions for a variety of problems including linear equations and matrix normal forms with the guiding design principle of re-usability. The most essential and driving design criterion for LinBox is that it is generic with respect to the domain of computation. This is because there are many and various representations of finite fields each of which is advantageous to use for some algorithm under some circumstance. The integral and rational number capabilities depend heavily on modular

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.

We consider the parallel iterative solution of indefinite linear systems given as augmented systems. Our numerical technique is based on an algebraic non-overlapping domain decomposition technique that only exploits the graph of the sparse matrix. This approach to high-performance, scalable solution of large sparse linear systems in parallel scientific computing, is to combine direct and iterative methods. We report numerical and parallel performance of the scheme on large matrices arising from the finite element discretization of linear elasticity in structural mechanics problems.

Parallel algebraic domain decomposition solver for the solution of augmented systems