This article presents a parallel, effective, and feature-complete recursive SPIKE algorithm that achieves near feature-parity with the standard linear algebra package banded linear system solver. First, we present a flexible parallel implementation of the recursive SPIKE scheme that aims at removing its original limitation that the number of cores/processors be restricted to powers of two. A new transpose solve option for SPIKE is then developed to satisfy a standard requirement of most numerical solver libraries. Finally, a pivoting recursive SPIKE strategy is presented as an alternative to the non-pivoting scheme to improve numerical stability. All these new enhancements lead to the release of a new black-box feature-complete SPIKE-OpenMP package that significantly improves upon the performance and scalability obtained with other state-of-the-art banded solvers.

We present efficient coding schemes and distributed implementations of erasure coded linear system solvers. Erasure coded computations belong to the class of algorithmic fault tolerance schemes. They are based on augmenting an input dataset, executing the algorithm on the augmented dataset, and in the event of a fault, recovering the solution from the corresponding augmented solution. This process can be viewed as the computational analog of erasure coded storage schemes. The proposed technique has a number of important benefits: (i) as the hardware platform scales in size and number of faults, our scheme yields increasing improvement in resource utilization, compared to traditional schemes; (ii) the proposed scheme is easy to code - the core algorithms remain the same; and (iii) the general scheme is flexible - accommodating a range of computation and communication tradeoffs. We present new coding schemes for augmenting the input matrix that satisfy the recovery equations of erasure coding with high probability in the event of random failures. These coding schemes also minimize fill (non-zero elements introduced by the coding block), while being amenable to efficient partitioning across processing nodes. We demonstrate experimentally that our scheme adds minimal overhead for fault tolerance, yields excellent parallel efficiency and scalability, and is robust to different fault arrival models.

Computational aeroacoustics (CAA) has emerged as a tool to complement theoretical and experimental approaches for robust and accurate prediction of sound levels from aircraft airframes and engines. CAA, unlike computational fluid dynamics (CFD), involves the accurate prediction of small-amplitude acoustic fluctuations and their correct propagation to the far field. In that respect, CAA poses significant challenges for researchers because the computational scheme should have high accuracy, good spectral resolution, and low dispersion and diffusion errors. A high-order compact finite difference scheme, which is implicit in space, can be used for such simulations because it fulfills the requirements for CAA. Usually, this method is parallelized using a transposition scheme; however, that approach has a high communication overhead. In this paper, we discuss the use of a parallel tridiagonal linear system solver based on the truncated SPIKE algorithm for reducing the communication overhead in our large eddy simulations. We report experimental results collected on two parallel computing platforms.

With availability of large-scale parallel platforms comprised of tens-of-thousands of processors and beyond, there is significant impetus for the development of scalable parallel sparse linear system solvers and preconditioners. An integral part of this design process, is the development of performance models capable of predicting performance and providing accurate cost models for the solvers and preconditioners. There has been some work in the past on characterizing performance of the iterative solvers themselves. In this paper, we investigate the problem of characterizing performance and scalability of banded preconditioners. Recent work has demonstrated the superior convergence properties and robustness of banded preconditioners, compared to state-of-the-art ILU family of preconditioners. Furthermore, when used in conjunction with efficient banded solvers, banded preconditioners are capable of significantly faster time-to solution. Our banded solver, the Truncated Spike algorithm is specifically designed for parallel performance and tolerance to deep memory hierarchies. Its regular structure is also highly amenable to accurate performance characterization. Using these characteristics, we derive the following results in this paper: (i) we develop parallel formulations of the Truncated Spike solver, (ii) we develop a highly accurate pseudo-analytical parallel performance model for our solver, (iii) we show excellent predication capabilities of our model - based on which we argue the high scalability of our solver. Our pseudo-analytical performance model is based on analytical performance characterization of each phase of our solver. These analytical models are then parameterized using actual runtime information on target platforms. An important consequence of our performance models is that they reveal underlying performance bottlenecks in both serial and parallel formulations. All of our results are validated on diverse heterogeneous multiclusters - platforms for which performance prediction is particularly challenging.

As the shared memory bus becomes a major performance bottleneck for many numerical applications on multicore chips, understanding how the increased parallelism on chip strains the memory bandwidth and hence affects the efficiency of parallel codes becomes a critical issue. This paper introduces the notion of memory access intensity to facilitate quantitative analysis of program's memory behavior on multicores which employ state-of-the-art prefetching hardware. Three numerical solvers for large scale sparse linear systems are used to demonstrate the estimation of memory access intensity and its effect on program performance.

Rapid evaluation of potentials in particle systems is an important and time-consuming step in many physical simulations. Over the past decade (1988-98), the development of treecodes such as the Fast Multipole Method (FMM) and the Barnes-Hut method has enabled large scale simulations in domains such as astrophysics, molecular dynamics, and material science. FMM and related methods rely on fixed degree polynomial (p) approximations of the potential of a set of points in a hierarchy. We present a sequence of results to illustrate that keeping the multipole degree constant can lead to large aggregate errors. An alternate strategy based on a careful selection of the multipole degree leads to asymptotically lower errors; while incurring minimal computation overhead for practical problem sizes. The paper presents theoretical results for computing the degree of a particle cluster interaction, the error associated with the interaction, the error associated with a particle for all of its interactions, and the computational complexity of the new method. These results show that it is possible to reduce the simulation error asymptotically while incurring minimal computational overhead. The paper also presents experimental validation of these results on a 32 processor Origin 2000 in the context of problems ranging from astrophysics to boundary element solvers. In addition to verifying theoretical results, we also show that it is possible to achieve excellent parallel speedup for the treecode.

Abstract: The problem of evaluating the potential due to a set of particles is an important and time- consuming one. The development of fast treecodes such as the Barnes-Hut and Fast Multipole Methods for n-body systems has enabled large scale simulations in astrophysics [9, 10, 13] and molecular dynamics [1]. Coupled with efficient parallel processing, these treecodes are capable of yielding several orders of magnitude improvement in performance [6, 14, 15]. In addition, treecodes have applications in the solution of dense linear systems arising from boundary element methods [3, 4, 5, 11, 12]. Using a p-term multipole expansion, the FMM reduces the complexity of a single timestep from O(n2) to O(p2n) and Barnes-Hut method reduces it to O(p2log n) for a uniform distribution. In this paper, we analyze the approximations introduced by these methods. We describe an algorithm that reduces the error significantly by selecting the multipole degree appropriately for different clusters. Furthermore, we show that for practical problem sizes, this increases the computational complexity marginally. We support our theoretical result with experiments in the context of particle simulations as well as boundary element methods. Our POSIX threads-based treecode yields excellent speedups on a 32 processor SGI Origin 2000, even for relatively small problems.

The method of moments is an important tool for solving boundary integral equations arising in a variety of applications. It transforms the physical problem into a dense linear system. Due to the large number of variables and the associated computational requirements, these systems are solved iteratively using methods such as GMRES, CG and its variants. The core operation of thes itertive solvers is the application of the system matrix to a vector. This requres O(n2) operations and memory using accurate dense methods. The computational complexity can be reduced to O(n log n) and the memory requirement to O(n) using hierarchical approximation techniques. The algorithmic speedup from approximation can be combined with parallelism to yield very fast dense solvers. In this paper, we present efficient parallel formulations of dense iterative solvers based on hierarchical approximations for solving the integral form of Laplace equation. We study the impact of various parameters on the accuracy and performance of the parallel solver. We present two preconditioning techniques for accelerating the convergence of the iterative solver. Thes techniques are based on an inner-outer scheme and a block diagonal scheme based on a truncated Green's function. We present detailed experimental results on up to 256 processors of a Cray T3D.

. Iterative schemes for solving sparse linear systems arising from elliptic PDEs are very suitable for efficient implementation on large scale multiprocessors. However, these methods rely heavily on effective preconditioners which must also be amenable to parallelization. In this paper, we present a novel method to obtain a preconditioned linear system which is solved using an iterative method. Each iteration comprises of a matrix-vector product with k sparse matrices (k log n), and can be computed in O(n) operations where n is the number of unknowns. The numerical convergence properties of our preconditioner are superior to the commonly used incomplete factorization preconditioners. Moreover, unlike the incomplete factorization preconditioners, our algorithm affords a higher degree of concurrency and doesn't require triangular system solves, thereby achieving the dual objective of good preconditioning and efficient parallel implementation. We describe our scheme for certain linear sy...

Matrix-vector products (mat-vecs) form the core of iterative methods used for solving dense linear systems. Often, these systems arise in the solution of integral equations used in electromagnetics, heat transfer, and wave propagation. In this paper, we present a parallel approximate method for computing mat-vecs used in the solution of integral equations. We use this method to compute dense mat-vecs of hundreds of thousands of elements. The combined speedups obtained from the use of approximate methods and parallel processing represent an improvement of several orders of magnitude over exact mat-vecs on uniprocessors. We demonstrate that our parallel formulation incurs minimal parallel processing overhead and scales up to a large number of processors. We study the impact of varying the accuracy of the approximate mat-vec on overall time and on parallel efficiency. Experimental results are presented for 256 processor Cray T3D and Thinking Machines CM5 parallel computers. We have achieved computation rates in excess of 5 GFLOPS on the T3D.

This paper analyzes the performance and scalability of an iteration of the preconditioned conjugate gradient algorithm on parallel architectures with a variety of interconnection networks, such as the mesh, the hypercube, and that of the CM-5 parallel computer. It is shown that for block-tridiagonal matrices resulting from two-dimensional finite difference grids, the communication overhead due to vector inner products dominates the communication overheads of the remainder of the computation on a large number of processors. However, with a suitable mapping, the parallel formulation of a PCG iteration is highly scalable for such matrices on a machine like the CM-5 whose fast control network practically eliminates the overheads due to inner product computation. The use of the truncated Incomplete Cholesky (IC) preconditioner can lead to further improvement in scalability on the CM-5 by a constant factor,as a result, a parallel formulation of the PCG algorithm with IC preconditioner may execute faster than that with a simple diagonal preconditioner even if the latter runs faster in a serial implementation.>

We present two new parallel formulations of the Barnes-Hut method. These parallel formulations are especially suited for simulations with irregular particle densities. We first present a parallel formulation that uses a static partitioning of the domain and assignment of subdomains to processors. We demonstrate that this scheme delivers acceptable load balance, and coupled with two collective communication operations, it yields good performance. We present a second parallel formulation which combines static decomposition of the domain with an assignment of subdomains to processors based on Morton ordering. This alleviates the load imbalance inherent in the first scheme. The second parallel formulation is inspired by two currently best known parallel algorithms for the Barnes-Hut method. We present an experimental evaluation of these schemes on a 256 processor nCUBE2 parallel computer for an astrophysical simulation.>

In this paper, we give an overview of the Cedar multiprocessor and present recent performance results. These include the performance of some computational kernels and the Perfect Benchmarks. We also present a methodology for judging parallel system performance and apply this methodology to Cedar, Cray YMP-8, and Thinking Machines CM-5.

The problem of finding an approximation of @@@@ = A†b (where A† is the pseudo-inverse of A ∈ @@@@m@@@@n with m ≥ n and rank(A) = n) is discussed. It is assumed that A is sparse but has neither a special pattern (as bandedness) nor a special property (as symmetry or positive definiteness). In this paper it is shown that preconditioners obtained by neglecting small elements during the decomposition of A into easily invertible matrices can be used efficiently with conjugate gradient-type methods if an adaptive strategy for deciding when an element is small is implemented. The resulting preconditioned methods are often better than the corresponding direct and pure iterative methods or those based on preconditionings in which elements are neglected when they appear in a predefined set of positions in the matrix, i.e. positional rather than numerical dropping. Numerical results are given to illustrate the performance of the CG-type methods preconditioned via numerical dropping.

An iterative method for the solution of nonsymmetric linear systems of equations is described and tested. The method, block symmetric successive over-relaxation with conjugate gradient acceleration (BSSOR), is remarkably robust and when applied to block tridiagonal systems allows parallelism in the computations. BSSOR compares favorably to unpreconditioned conjugate gradient-like algorithms in efficiency, and although generally slower than preconditioned methods it is far more reliable. The concept behind BSSOR can, in general, be applied to sparse linear systems (even if they are singular), sparse nonlinear systems of equations and least squares problems.

No abstract available.

We present an algorithm for solving banded positive defimte linear systems on a multiprocessor computer whose number of processors p is much less than the order of the system n.Assuming that the banded matrix, of bandwidth 2m + 1, is stored in the global memory by diagonals as several onedlmensmnal arrays, we consider the time required by several alignment networks for allocating the appropriate data to the local memory of each processor.We demonstrate that the time required in this preprocessmg stage does not exceed that required by the algorithm provided we use a shuffle exchange, a plpelmed shuffle exchange, or a crossbar switch.Once the data are allocated in the local memories, the algorithm requires only a "nearest neighbor" alignment network to achieve a total time of O(m'~n/p).The total cost of the algorithm is minimized when p ~ ~.

In this paper three stable parallel algorithms for solving dense and tndlagonai systems of lmear equations are discussed The algorithms are based on Givens' reduction of a matrix to the upper triangular form The algorithm for the dense case requires O(n) time steps compared to O(n log n) steps for Gausslan ehmmatlon with pivoting (in the absence of certain features of machine logic and hardware) For the trldlagonal case, one of the algorithms presented here is superior to the best previous algorithm in that with a modest increase in time It does not fall if any of the leading pnnclpal submatrlces is singular, the probablhty of over-or underflow is minimized, and the error bound does not grow exponentially Furthermore, it is most statable when only a hmtted number of processors ts available.

article Free Access Share on Practical Parallel Band Triangular System Solvers Authors: S. C. Chen Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, IL Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, ILView Profile , D. J. Kuck Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, IL Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, ILView Profile , A. H. Sameh Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, IL Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, ILView Profile Authors Info & Claims ACM Transactions on Mathematical SoftwareVolume 4Issue 3Sept. 1978 pp 270–277https://doi.org/10.1145/355791.355797Published:01 September 1978Publication History 80citation447DownloadsMetricsTotal Citations80Total Downloads447Last 12 Months20Last 6 weeks3 Get Citation AlertsNew Citation Alert added!This alert has been successfully added and will be sent to:You will be notified whenever a record that you have chosen has been cited.To manage your alert preferences, click on the button below.Manage my Alerts New Citation Alert!Please log in to your account Save to BinderSave to BinderCreate a New BinderNameCancelCreateExport CitationPublisher SiteeReaderPDF

EffectsA study of relative error propagation In numerical algorithms is given.Roundoff errors and their effects are Identified, and a method IS presented which efficiently calculates these effects Present methods in use are O(l 2) The time and storage required by the method presented here is O(l), where l Is the number of operations in the straight-line algorithm being analyzed These savings will enable the analysis of larger and more complex algorithms Key Words and Phrases automatic roundoff analysis, numerical stability, numerical linear algebra CR Categories: 5.10, 5 11, 5 14

The error properties of floating-point arithmetic using various rounding methods (including ROM rounding, a new scheme) are analyzed. Guard digits are explained, and the rounding schemes' effectiveness are evaluated and compared.

We show that if the size of the tridiagonal matrix in any given iteration is n, then the parallel QR algorithm requires 0(log2n) steps with 0(n) processors per iteration and no square roots. This results in a speedup of 0(n/log2n) over the sequential algorithm with an efficiency of 0(1/log2n). We also give an error analysis of the parallel triangular system solvers used in each iteration.

ROM-rounding is introduced and is shown to compare favorably with existing floating-point rounding methods on design considerations and on performance over a series of error tests. The error-retarding value of guard digits, of rounding the aligned operand, and of rounding in general are discussed.

Input/output economic models lead to large dense unsymmetric systems of linear equations. While solutions of the equations are of primary interest, the situation is somewhat unusual in that the elements of the Leontief inverse are physically meaningful and hence of value. In the analysis of such systems we would like to determine the effects of the solution of the Leontief equations due to changes in the transaction data. In particular we discuss modifications down a column of the direct coefficients matrix, modifications along a row, and disaggregation of the industries in the new system. Numerical technique are available in the literature for obtaining new solutions due to these changes with far less operations than those required to solve the system ‘from scratch’.