A tearing-based hybrid parallel banded linear system solver

A new parallel algorithm for the solution of banded linear systems is proposed. The scheme tears the coefficient matrix into several overlapped independent blocks in which the size of the overlap is equal to the system’s bandwidth. A corresponding splitting of the right-hand side is also provided. The resulting independent, and smaller size, linear systems are solved under the constraint that the solutions corresponding to the overlap regions are identical. This results in a linear system whose size is proportional to the sum of the overlap regions which we refer to as the “balance” system. We propose a solution strategy that does not require obtaining this “balance” system explicitly. Once the balance system is solved, retrieving the rest of the solution can be realized with almost perfect parallelism. Our proposed algorithm is a hybrid scheme that combines direct and iterative methods for solving a single banded system of linear equations on parallel architectures. It has broad applications in finite-element analysis, particularly as a parallel solver of banded preconditioners that can be used in conjunction with outer Krylov iterative schemes.     

Solving sparse linear systems by an abs-method that corresponds toLU-decomposition

Abaffy, Broyden, and Spedicato (ABS) have recently proposed a class of direct methods for solving nonsparse linear systems. It is the purpose of this paper to demonstrate that with proper choice of parameters, ABS methods exploit sparsity in a natural way. In particular, we study the choice of parameters which corresponds to an LU-decomposition of the coefficient matrix. In many cases, the fill-in, represented by the nonzero elements of the deflection matrix, uses less storage than the corresponding fill-in of an explicit LU factorization. Hence the above can be a viable and simple method for solving sparse linear systems. A simple reordering scheme for the coefficient matrix is also proposed for the purpose of reducing fill-in of the deflection matrices.     

on hybrid preconditioning methods for large sparse saddle-point problems

Based on the block-triangular product approximation to a 2-by-2 block matrix, a class of hybrid preconditioning methods is designed for accelerating the MINRES method for solving saddle-point problems. The appropriate values for the parameters involved in the new preconditioners are estimated, so that the numerical conditioning and the spectral property of the saddle-point matrix of the linear system can be substantially improved. Several practical hybrid preconditioners and the corresponding preconditioning iterative methods are constructed and studied, too.     

Wavelet sparse approximate inverse preconditioners

We show how to use wavelet compression ideas to improve the performance of approximate inverse preconditioners. Our main idea is to first transform the inverse of the coefficient matrix into a wavelet basis, before applying standard approximate inverse techniques. In this process, smoothness in the entries ofA ⁻¹ are converted into small wavelet coefficients, thus allowing a more efficient approximate inverse approximation. We shall justify theoretically and numerically that our approach is effective for matrices with smooth inverses. Supported by grants from ONR: ONR-N00014-92-J-1890, and the Army Research Office: DAAL-03-91-C-0047 (Univ. of Tenn. subcontract ORA4466.04 Amendment 1 and 2). The first and the third author also acknowledge support from RIACS/NASA Ames NAS 2-96027 and the Alfred P. Sloan Foundation as Doctoral Dissertation Fellows, respectively. the work was supported by the Natural Sciences and Engineering Research Council of Canada, the Information Technology Research Centre (which is funded by the Province of Ontario), and RIACS/NASA Ames NAS 2-96027.     

A domain-decomposing parallel sparse linear system solver

The solution of large sparse linear systems is often the most time-consuming part of many science and engineering applications. Computational fluid dynamics, circuit simulation, power network analysis, and material science are just a few examples of the application areas in which large sparse linear systems need to be solved effectively. In this paper, we introduce a new parallel hybrid sparse linear system solver for distributed memory architectures that contains both direct and iterative components. We show that by using our solver one can alleviate the drawbacks of direct and iterative solvers, achieving better scalability than with direct solvers and more robustness than with classical preconditioned iterative solvers. Comparisons to well-known direct and iterative solvers on a parallel architecture are provided.     

A domain-decomposing parallel sparse linear system solver

The solution of large sparse linear systems is often the most time-consuming part of many science and engineering applications. Computational fluid dynamics, circuit simulation, power network analysis, and material science are just a few examples of the application areas in which large sparse linear systems need to be solved effectively. In this paper, we introduce a new parallel hybrid sparse linear system solver for distributed memory architectures that contains both direct and iterative components. We show that by using our solver one can alleviate the drawbacks of direct and iterative solvers, achieving better scalability than with direct solvers and more robustness than with classical preconditioned iterative solvers. Comparisons to well-known direct and iterative solvers on a parallel architecture are provided.     

On structure-oriented hybrid two-stage iteration methods for the large and sparse blocked system of linear equations

In this paper, we first present a class of structure-oriented hybrid two-stage iteration methods for solving the large and sparse blocked system of linear equations, as well as the saddle point problem as a special case. And the new methods converge to the solution under suitable restrictions, for instance, when the coefficient matrix is positive stable matrix generally. Numerical experiments for a model generalized saddle point problem are given, and the results show that our new methods are feasible and efficient, and converge faster than the Classical Uzawa Method.     

On CSCS-based iteration methods for Toeplitz system of weakly nonlinear equations

For Toeplitz system of weakly nonlinear equations, by using the separability and strong dominance between the linear and the nonlinear terms and using the circulant and skew-circulant splitting (CSCS) iteration technique, we establish two nonlinear composite iteration schemes, called Picard-CSCS and nonlinear CSCS-like iteration methods, respectively. The advantage of these methods is that they do not require accurate computation and storage of Jacobian matrix, and only need to solve linear sub-systems of constant coefficient matrices. Therefore, computational workloads and computer storage may be saved in actual implementations. Theoretical analysis shows that these new iteration methods are local convergent under suitable conditions. Numerical results show that both Picard-CSCS and nonlinear CSCS-like iteration methods are feasible and effective for some cases.     

A parallel implementation of the restarted GMRES iterative algorithm for nonsymmetric systems of linear equations

We describe the parallelisation of the GMRES(c) algorithm and its implementation on distributed-memory architectures, using both networks of transputers and networks of workstations under the PVM message-passing system. The test systems of linear equations considered are those derived from five-point finite-difference discretisations of partial differential equations. A theoretical model of the computation and communication phases is presented which allows us to decide for which values of the parameterc our implementation executes efficiently. The results show that for reasonably large discretisation grids the implementations are effective on a large number of processors.     

A preconditioning approach to the pagerank computation problem

Some spectral properties of the transition matrix of an oriented graph indicate the preconditioning of Euler-Richardson (ER) iterative scheme as a good way to compute efficiently the vertexrank vector associated with such graph. We choose the preconditioner from an algebra U of matrices, thereby obtaining an ER_U method, and we observe that ER_U can outperform ER in terms of rate of convergence. The proposed preconditioner can be updated at a very low cost whenever the graph changes, as is the case when it represents a generic set of information. The particular U utilized requires a surplus of operations per step and memory allocations, which make ER_U superior to ER for not too wide graphs. However, the observed high improvement in convergence rate obtained by preconditioning and the general theory developed, are a reason for investigating different choices of U, more efficient for huge graphs.     

A hybrid Arnoldi-Faber iterative method for nonsymmetric systems of linear equations

Summary We present here a new hybrid method for the iterative solution of large sparse nonsymmetric systems of linear equations, say of the formAx=b, whereA ^{N, N} , withA nonsingular, andb ^N are given. This hybrid method begins with a limited number of steps of the Arnoldi method to obtain some information on the location of the spectrum ofA, and then switches to a Richardson iterative method based on Faber polynomials. For a polygonal domain, the Faber polynomials can be constructed recursively from the parameters in the Schwarz-Christoffel mapping function. In four specific numerical examples of non-normal matrices, we show that this hybrid algorithm converges quite well and is approximately as fast or faster than the hybrid GMRES or restarted versions of the GMRES algorithm. It is, however, sensitive (as other hybrid methods also are) to the amount of information on the spectrum ofA acquired during the first (Arnoldi) phase of this procedure.     

Estimating the largest singular values of large sparse matrices via modified moments

We describe a procedure for determining a few of the largest singular values of a large sparse matrix. The method by Golub and Kent which uses the method of modified moments for estimating the eigenvalues of operators used in iterative methods for the solution of linear systems of equations is appropriately modified in order to generate a sequence of bidiagonal matrices whose singular values approximate those of the original sparse matrix. A simple Lanczos recursion is proposed for determining the corresponding left and right singular vectors. The potential asynchronous computation of the bidiagonal matrices using modified moments with the iterations of an adapted Chebyshev semi-iterative (CSI) method is an attractive feature for parallel computers. Comparisons in efficiency and accuracy with an appropriate Lanczos algorithm (with selective re-orthogonalization) are presented on large sparse (rectangular) matrices arising from applications such as information retrieval and seismic reflection tomography. This procedure is essentially motivated by the theory of moments and Gauss quadrature.This author's work was supported by the National Science Foundation under grants NSF CCR-8717492 and CCR-910000N (NCSA), the U.S. Department of Energy under grant DOE DE-FG02-85ER25001, and the Air Force Office of Scientific Research under grant AFOSR-90-0044 while at the University of Illinois at Urbana-Champaign Center for Supercomputing Research and Development.This author's work was supported by the U.S. Army Research Office under grant DAAL03-90-G-0105, and the National Science Foundation under grant NSF DCR-8412314.     

Computing projections with LSQR

LSQR uses the Golub-Kahan bidiagonalization process to solve sparse least-squares problems with and without regularization. In some cases, projections of the right-hand side vector are required, rather than the least-squares solution itself. We show that projections may be obtained from the bidiagonalization as linear combinations of (the-oretically) orthogonal vectors. Even the least-squares solution may be obtained from orthogonal vectors, perhaps more accurately than the usual LSQR solution. (However, LSQR has proved equally good in all examples so far.) Presented at the Cornelius Lanczos International Centenary Conference, North Carolina State University, Raleigh, NC, December 1993. Partially supported by Department of Energy grant DE-FG03-92ER25117, National Science Foundation grants DMI-9204208 and DMI-9500668, and Office of Naval Research grants N00014-90-J-1242 and N00014-96-1-0274.     

Two-dimensional block partitionings for the parallel sparse Cholesky factorization

This paper presents a discussion on 2D block mappings for the sparse Cholesky factorization on parallel MIMD architectures with distributed memory. It introduces the fan-in algorithm in a general manner and proposes several mapping strategies. The grid mapping with row balancing, inspired by Rothberg's work (1994), is proved to be more robust than the original fan-out algorithm. Even more efficient is the proportional mapping, as shown by the experiments on a 32 processor IBM SP1 and on a Cray T3D. Subforest-to-subcube mappings are also considered and give good results on the T3D. This revised version was published online in August 2006 with corrections to the Cover Date.     

Using domain decomposition to find graph bisectors

In this paper we introduce a three-step approach to find a vertex bisector of a graph. The first step finds adomain decomposition of the graph, consisting of a set of domains and a multisector. Eachdomain is a connected subgraph, and themultisector contains the remaining vertices that separate the domains from each other. The second step uses a block variant of the Kernighan-Lin scheme to find a bisector that is a subset of the multisector. The third step improves the bisector by bipartite graph matching. Experimental results show this domain decomposition method finds graph partitions that compare favorably with a state-of-the-art multilevel partitioning scheme in both quality and execution time. This research was supported in part by the ARPA Contract DABT63-95-C-0122, and in part by the Natural Sciences and Engineering Research Council of Canada under grant A5509.     

A numerical study of optimized sparse preconditioners

Preconditioning strategies based on incomplete factorizations and polynomial approximations are studied through extensive numerical experiments. We are concerned with the question of the optimal rate of convergence that can be achieved for these classes of preconditioners.Our conclusion is that the well-known Modified Incomplete Cholesky factorization (MIC), cf. e.g., Gustafsson [20], and the polynomial preconditioning based on the Chebyshev polynomials, cf. Johnson, Micchelli and Paul [22], have optimal order of convergence as applied to matrix systems derived by discretization of the Poisson equation. Thus for the discrete two-dimensional Poisson equation withn unknowns,O(n ^1/4) andO(n ^1/2) seem to be the optimal rates of convergence for the Conjugate Gradient (CG) method using incomplete factorizations and polynomial preconditioners, respectively. The results obtained for polynomial preconditioners are in agreement with the basic theory of CG, which implies that such preconditioners can not lead to improvement of the asymptotic convergence rate.By optimizing the preconditioners with respect to certain criteria, we observe a reduction of the number of CG iterations, but the rates of convergence remain unchanged.Supported by The Norwegian Research Council for Science and the Humanities (NAVF) under grants no. 413.90/002 and 412.93/005.Supported by The Royal Norwegian Council for Scientific and Industrial Research (NTNF) through program no. STP.28402: Toolkits in industrial mathematics.     

A fast solver for the first biharmonic boundary value problem

Summary This paper provides a fast and storage-saving method for the solution of the first biharmonic boundary value problem (b.v.p.). The b.v.p. is approximated via a special variational finite difference technique suggested earlier by V.G. Korneev. It is shown theoretically that our method produces an approximate solution to the finite difference equations inO(NlnNln^–1) arithmetical operations, whereN is the number of unknowns and (0<<1) denotes the relative accuracy required. The numerical results obtained by our computer code CGMFC decisively substantiate the theoretical estimates given.     

A class of Uzawa-SOR methods for saddle point problems

In this paper, we consider a class of Uzawa-SOR methods for saddle point problems, and prove the convergence of the proposed methods. We solve a lower triangular system per iteration in the proposed methods, instead of solving a linear equation Az=b. Actually, the new methods can be considered as an inexact iteration method with the Uzawa as the outer iteration and the SOR as the inner iteration. Although the proposed methods cannot achieve the same convergence rate as the GSOR methods proposed by Bai et al. [Z.-Z. Bai, B.N. Parlett, Z.-Q. Wang, On generalized successive overrelaxation methods for augmented linear systems, Numer. Math. 102 (2005) 1-38], but our proposed methods have less workloads per iteration step. Experimental results show that our proposed methods are feasible and effective.     

A practical termination criterion for the conjugate gradient method

The conjugate gradient method for the iterative solution of a set of linear equationsAx=b is essentially equivalent to the Lanczos method, which implies that approximations to certain eigen-values ofA can be obtained at low cost. In this paper it is shown how the smallest active eigenvalue ofA can be cheaply approximated, and the usefulness of this approximation for a practical termination criterion for the conjugate gradient method is studied. It is proved that this termination criterion is reliable in many relevant situations.