A tearing-based hybrid parallel sparse linear system solver

We propose a hybrid sparse system solver for handling linear systems using algebraic domain decomposition-based techniques. The solver consists of several stages. The first stage uses a reordering scheme that brings as many of the largest matrix elements as possible closest to the main diagonal. This is followed by partitioning the coefficient matrix into a set of overlapped diagonal blocks that contain most of the largest elements of the coefficient matrix. The only constraint here is to minimize the size of each overlap. Separating these blocks into independent linear systems with the constraint of matching the solution parts of neighboring blocks that correspond to the overlaps, we obtain a balance system. This balance system is not formed explicitly and has a size that is much smaller than the original system. Our novel solver requires only a one-time factorization of each diagonal block, and in each outer iteration, obtaining only the upper and lower tips of a solution vector where the size of each tip is equal to that of the individual overlap. This scheme proves to be scalable on clusters of nodes in which each node has a multicore architecture. Numerical experiments comparing the scalability of our solver with direct and preconditioned iterative methods are also presented.     

An incomplete-factorization preconditioning using repeated red-black ordering

Summary The convergence of the conjugate gradient method for the iterative solution of large systems of linear equations depends on proper preconditioning matrices. We present an efficient incomplete-factorization preconditioning based on a specific, repeated red-black ordering scheme and cyclic reduction. For the Dirichlet model problem, we prove that the condition number increases asymptotically slower with the number of equations than for usual incomplete factorization methods. Numerical results for symmetric and non-symmetric test problems and on locally refined grids demonstrate the performance of this method, especially for large linear systems.     

On the convergence of splittings for semidefinite linear systems

Recently, Lee et al. [Young-ju Lee, Jinbiao Wu, Jinchao Xu, Ludmil Zikatanov, On the convergence of iterative methods for semidefinite linear systems, SIAM J. Matrix Anal. Appl. 28 (2006) 634-641] introduce new criteria for the semi-convergence of general iterative methods for semidefinite linear systems based on matrix splitting. The new conditions generalize the classical notion of P-regularity introduced by Keller [H.B. Keller, On the solution of singular and semidefinite linear systems by iterations, SIAM J. Numer. Anal. 2 (1965) 281-290]. In view of their results, we consider here stipulations on a splitting A=M-N, which lead to fixed point systems such that, the iterative scheme converges to a weighted Moore-Penrose solution to the system Ax=b. Our results extend the result of Lee et al. to a more general case and we also show that it requires less restrictions on the splittings than Keller’s P-regularity condition to ensure the convergence of iterative scheme.     

A stationary iterative pseudoinverse algorithm

Iterative methods applied to the normal equationsA ^T Ax=A ^T b are sometimes used for solving large sparse linear least squares problems. However, when the matrix is rank-deficient many methods, although convergent, fail to produce the unique solution of minimal Euclidean norm. Examples of such methods are the Jacobi and SOR methods as well as the preconditioned conjugate gradient algorithm. We analyze here an iterative scheme that overcomes this difficulty for the case of stationary iterative methods. The scheme combines two stationary iterative methods. The first method produces any least squares solution whereas the second produces the minimum norm solution to a consistent system. This work was supported by the Swedish Research Council for Engineering Sciences, TFR.     

A generalized successive overrelaxation method for least squares problems

In this paper a new iterative method is given for solving large sparse least squares problems and computing the minimum norm solution to underdetermined consistent linear systems. The new scheme is called the generalized successive overrelaxation (GSOR) method and is shown to be convergent ifA is full column rank. The GSOR method involves a parameter ρ and an auxiliary matrixP. One can choose matrix P so that the GSOR method only involves matrix and vector operations; therefore the GSOR method is suitable for parallel computations. Besides, the GSOR method can be combined with preconditioning techniques, and therefore can be expected to be more effective. This author's work was supported by Natural Science Foundation of Liaoning Province, China.     

Schwarz alternating and iterative refinement methods for mixed formulations of elliptic problems,part I: Algorithms and numerical results

Summary We describe sequential and parallel algorithms based on the Schwarz alternating method for the solution of mixed finite element discretizations of elliptic problems using the Raviart-Thomas finite element spaces. These lead to symmetric indefinite linear systems and the algorithms have some similarities with the traditional block Gauss-Seidel or block Jacobi methods with overlapping blocks. The indefiniteness requires special treatment. The sub-blocks used in the algorithm correspond to problems on a coarse grid and some overlapping subdomains and is based on a similar partition used in an algorithm of Dryja and Widlund for standard elliptic problems. If there is sufficient overlap between the subdomains, the algorithm converges with a rate independent of the mesh size, the number of subdomains and discontinuities of the coefficients. Extensions of the above algorithms to the case of local grid refinement is also described. Convergence theory for these algorithms will be presented in a subsequent paper.This work was supported in part by the National Science Foundation under Grant NSF-CCR-8903003, while the author was a graduate student at New York University, and in part by the Army Research Office under Grant DAAL 03-91-G-0150, while the author was a Visiting Assistant Researcher at UCLA     

Minimum residual methods for augmented systems

For large systems of linear equations, iterative methods provide attractive solution techniques. We describe the applicability and convergence of iterative methods of Krylov subspace type for an important class of symmetric and indefinite matrix problems, namely augmented (or KKT) systems. Specifically, we consider preconditioned minimum residual methods and discuss indefinite versus positive definite preconditioning. For a natural choice of starting vector we prove that when the definite and indenfinite preconditioners are related in the obvious way, MINRES (which is applicable in the case of positive definite preconditioning) and full GMRES (which is applicable in the case of indefinite preconditioning) give residual vectors with identical Euclidean norm at each iteration. Moreover, we show that the convergence of both methods is related to a system of normal equations for which the LSQR algorithm can be employed. As a side result, we give a rare example of a non-trivial normal(1) matrix where the corresponding inner product is explicitly known: a conjugate gradient method therefore exists and can be employed in this case. This work was supported by British Council/German Academic Exchange Service Research Collaboration Project 465 and NATO Collaborative Research Grant CRG 960782     

Application of a block modified Chebyshev algorithm to the iterative solution of symmetric linear systems with multiple right hand side vectors

Summary. An adaptive Richardson iteration method is described for the solution of large sparse symmetric positive definite linear systems of equations with multiple right-hand side vectors. This scheme ``learns' about the linear system to be solved by computing inner products of residual matrices during the iterations. These inner products are interpreted as block modified moments. A block version of the modified Chebyshev algorithm is presented which yields a block tridiagonal matrix from the block modified moments and the recursion coefficients of the residual polynomials. The eigenvalues of this block tridiagonal matrix define an interval, which determines the choice of relaxation parameters for Richardson iteration. Only minor modifications are necessary in order to obtain a scheme for the solution of symmetric indefinite linear systems with multiple right-hand side vectors. We outline the changes required. Received April 22, 1993     

A parallel balance scheme for banded linear systems

A parallel algorithm is proposed for the solution of narrow banded non‐symmetric linear systems. The linear system is partitioned into blocks of rows with a small number of unknowns common to multiple blocks. Our technique yields a reduced system defined only on these common unknowns which can then be solved by a direct or iterative method. A projection based extension to this approach is also proposed for computing the reduced system implicitly, which gives rise to an inner–outer iteration method. In addition, the product of a vector with the reduced system matrix can be computed efficiently on a multiprocessor by concurrent projections onto subspaces of block rows. Scalable implementations of the algorithm can be devized for hierarchical parallel architectures by exploiting the two‐level parallelism inherent in the method. Our experiments indicate that the proposed algorithm is a robust and competitive alternative to existing methods, particularly for difficult problems with strong indefinite symmetric part. Copyright © 2001 John Wiley & Sons, Ltd.     

Fast inexact subspace iteration for generalized eigenvalue problems with spectral transformation

We study inexact subspace iteration for solving generalized non-Hermitian eigenvalue problems with spectral transformation, with focus on a few strategies that help accelerate preconditioned iterative solution of the linear systems of equations arising in this context. We provide new insights into a special type of preconditioner with “tuning” that has been studied for this algorithm applied to standard eigenvalue problems. Specifically, we propose an alternative way to use the tuned preconditioner to achieve similar performance for generalized problems, and we show that these performance improvements can also be obtained by solving an inexpensive least squares problem. In addition, we show that the cost of iterative solution of the linear systems can be further reduced by using deflation of converged Schur vectors, special starting vectors constructed from previously solved linear systems, and iterative linear solvers with subspace recycling. The effectiveness of these techniques is demonstrated by numerical experiments.     

Optimum first and second order extrapolations of Successive Overrelaxation type methods for certain classes of matrices

This paper deals with the iterative solution of the linear systemx=Bx+c when its Jacobi matrixB is weakly 2-cyclic consistently ordered and has a complex eigenvalue spectrum which lies on a straight-line segment. The optimization problem of the following three methods is considered and solved: i) The extrapolation of the optimum Successive Overrelaxation (SOR) ii) The second order extrapolation of a good SOR and iii) The second order extrapolation of the Gauss-Seidel method. In addition a variant of the second order methods considered, suitable for the solution of the system even ifB isnot necessarily weakly 2-cyclic consistently ordered, is proposed. Finally a reference to a theoretical comparison of the various optimum methods in the paper is made and their asymptotic convergence factors for selected eigenvalue spectra are illustrated in a Table in support of the theory developed.     

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.     

Augmented conjugate gradient. Application in an iterative process for the solution of scattering problems

We discuss the application of an augmented conjugate gradient to the solution of a sequence of linear systems of the same matrix appearing in an iterative process for the solution of scattering problems. The conjugate gradient method applied to the first system generates a Krylov subspace, then for the following systems, a modified conjugate gradient is applied using orthogonal projections on this subspace to compute an initial guess and modified descent directions leading to a better convergence. The scattering problem is treated via an Exact Controllability formulation and a preconditioned conjugate gradient algorithm is introduced. The set of linear systems to be solved are associated to this preconditioning. The efficiency of the method is tested on different 3D acoustic problems. This revised version was published online in August 2006 with corrections to the Cover Date.     

A communication-less parallel algorithm for tridiagonal Toeplitz systems

Diagonally dominant tridiagonal Toeplitz systems of linear equations arise in many application areas and have been well studied in the past. Modern interest in numerical linear algebra is often focusing on solving classic problems in parallel. In McNally [Fast parallel algorithms for tri-diagonal symmetric Toeplitz systems, MCS Thesis, University of New Brunswick, Saint John, 1999], an m processor Split & Correct algorithm was presented for approximating the solution to a symmetric tridiagonal Toeplitz linear system of equations. Nemani [Perturbation methods for circulant-banded systems and their parallel implementation, Ph.D. Thesis, University of New Brunswick, Saint John, 2001] and McNally (2003) adapted the works of Rojo [A new method for solving symmetric circulant tri-diagonal system of linear equations, Comput. Math. Appl. 20 (1990) 61–67], Yan and Chung [A fast algorithm for solving special tri-diagonal systems, Computing 52 (1994) 203–211] and McNally et al. [A split-correct parallel algorithm for solving tri-diagonal symmetric Toeplitz systems, Internat. J. Comput. Math. 75 (2000) 303–313] to the non-symmetric case. In this paper we present relevant background from these methods and then introduce an m processor scalable communication-less approximation algorithm for solving a diagonally dominant tridiagonal Toeplitz system of linear equations.     

A cascadic multigrid algorithm for semilinear elliptic problems

Summary. We propose a cascadic multigrid algorithm for a semilinear elliptic problem. The nonlinear equations arising from linear finite element discretizations are solved by Newton's method. Given an approximate solution on the coarsest grid on each finer grid we perform exactly one Newton step taking the approximate solution from the previous grid as initial guess. The Newton systems are solved iteratively by an appropriate smoothing method. We prove that the algorithm yields an approximate solution within the discretization error on the finest grid provided that the start approximation is sufficiently accurate and that the initial grid size is sufficiently small. Moreover, we show that the method has multigrid complexity. Received February 12, 1998 / Revised version received July 22, 1999 / Published online June 8, 2000     

Iterative refinement enhances the stability ofQR factorization methods for solving linear equations

Iterative refinement is a well-known technique for improving the quality of an approximate solution to a linear system. In the traditional usage residuals are computed in extended precision, but more recent work has shown that fixed precision is sufficient to yield benefits for stability. We extend existing results to show that fixed precision iterative refinement renders anarbitrary linear equations solver backward stable in a strong, componentwise sense, under suitable assumptions. Two particular applications involving theQR factorization are discussed in detail: solution of square linear systems and solution of least squares problems. In the former case we show that one step of iterative refinement suffices to produce a small componentwise relative backward error. Our results are weaker for the least squares problem, but again we find that iterative refinement improves a componentwise measure of backward stability. In particular, iterative refinement mitigates the effect of poor row scaling of the coefficient matrix, and so provides an alternative to the use of row interchanges in the HouseholderQR factorization. A further application of the results is described to fast methods for solving Vandermonde-like systems.     

Time-harmonic solution for acousto-elastic interaction with controllability and spectral elements

The classical way of solving the time-harmonic linear acousto-elastic wave problem is to discretize the equations with finite elements or finite differences. This approach leads to large-scale indefinite complex-valued linear systems. For these kinds of systems, it is difficult to construct efficient iterative solution methods. That is why we use an alternative approach and solve the time-harmonic problem by controlling the solution of the corresponding time dependent wave equation.In this paper, we use an unsymmetric formulation, where fluid-structure interaction is modeled as a coupling between pressure and displacement. The coupled problem is discretized in space domain with spectral elements and in time domain with central finite differences. After discretization, exact controllability problem is reformulated as a least-squares problem, which is solved by the conjugate gradient method.     

Parallel factorization of banded linear matrices using a systolic array processor

A parallel algorithm for calculating theQR factorization of a banded system of linear equations using a systolic array processor is presented and an application to spline fitting is given. The major advantage of the method is that the size of the processor array is fixed by the size of the bandwidth of the system to be solved. This allows the factorization of large systems using small systolic arrays. The cost of the method, in terms of storage and time, is optimal.     

AIR Tools — A MATLAB package of algebraic iterative reconstruction methods

We present a MATLAB package with implementations of several algebraic iterative reconstruction methods for discretizations of inverse problems. These so-called row action methods rely on semi-convergence for achieving the necessary regularization of the problem. Two classes of methods are implemented: Algebraic Reconstruction Techniques (ART) and Simultaneous Iterative Reconstruction Techniques (SIRT). In addition we provide a few simplified test problems from medical and seismic tomography. For each iterative method, a number of strategies are available for choosing the relaxation parameter and the stopping rule. The relaxation parameter can be fixed, or chosen adaptively in each iteration; in the former case we provide a new “training” algorithm that finds the optimal parameter for a given test problem. The stopping rules provided are the discrepancy principle, the monotone error rule, and the NCP criterion; for the first two methods “training” can be used to find the optimal discrepancy parameter.     

Incomplete factorizations of singular linear systems

Recent works have shown that, whenA is a Stieltjes matrix, its so-called modified incomplete factorizations provide effective preconditioning matrices for solvingAx=b by polynomially accelerated iterative methods. We extend here these results to the singular case with the conclusion that the latter techniques are able to solve singular systems at the same speed as regular systems.Research supported by the Fonds National de la Recherche Scientifique (Belgium) — Aspirant.