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.     

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.     

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.     

On deflation and singular symmetric positive semi-definite matrices

For various applications, it is well-known that the deflated ICCG is an efficient method for solving linear systems with invertible coefficient matrix. We propose two equivalent variants of this deflated ICCG which can also solve linear systems with singular coefficient matrix, arising from discretization of the discontinuous Poisson equation with Neumann boundary conditions. It is demonstrated both theoretically and numerically that the resulting methods accelerate the convergence of the iterative process.     

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.     

Factoring matrices with a tree-structured sparsity pattern

Let A be a matrix whose sparsity pattern is a tree with maximal degree d_max. We show that if the columns of A are ordered using minimum degree on |A|+|A|^∗, then factoring A using a sparse LU with partial pivoting algorithm generates only O(d_maxn) fill, requires only O(d_maxn) operations, and is much more stable than LU with partial pivoting on a general matrix. We also propose an even more efficient and just-as-stable algorithm called sibling-dominant pivoting. This algorithm is a strict partial pivoting algorithm that modifies the column preordering locally to minimize fill and work. It leads to only O(n) work and fill. More conventional column pre-ordering methods that are based (usually implicitly) on the sparsity pattern of |A|^∗|A| are not as efficient as the approaches that we propose in this paper.     

Computing permanents via determinants for some classes of sparse matrices

Starting from recent formulas for calculating the permanents of some sparse circulant matrices, we obtain more general formulas expressing the permanents of a wider class of matrices as a linear combination of appropriate determinants.     

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.     

Graph properties for splitting with grounded Laplacian matrices

With any undirected, connected graphG containing no self-loops one can associate the Laplacian matrixL(G). It is also the (singular) admittance matrix of a resistive network with all conductances taken to be unity. While solving the linear system involved, one of the vertices is grounded, so the coefficient matrix is a principal submatrix ofL which we will call the grounded Laplacian matrixL ₁. In this paper we consider iterative solutions of such linear systems using certain regular splittings ofL ₁ and derive an upper bound for the spectral radius of the iteration matrix in terms of the properties of the graphG.This work was supported by the Academy of Finland     

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.     

On generalized symmetric SOR method for augmented systems

In this paper, we establish the generalized symmetric SOR method (GSSOR) for solving the large sparse augmented systems of linear equations, which is the extension of the SSOR iteration method. The convergence of the GSSOR method for augmented systems is studied. Numerical resume shows that this method is effective.     

Approximate iterations for structured matrices

Important matrix-valued functions f (A) are, e.g., the inverse A ⁻¹, the square root and the sign function. Their evaluation for large matrices arising from pdes is not an easy task and needs techniques exploiting appropriate structures of the matrices A and f (A) (often f (A) possesses this structure only approximately). However, intermediate matrices arising during the evaluation may lose the structure of the initial matrix. This would make the computations inefficient and even infeasible. However, the main result of this paper is that an iterative fixed-point like process for the evaluation of f (A) can be transformed, under certain general assumptions, into another process which preserves the convergence rate and benefits from the underlying structure. It is shown how this result applies to matrices in a tensor format with a bounded tensor rank and to the structure of the hierarchical matrix technique. We demonstrate our results by verifying all requirements in the case of the iterative computation of A ⁻¹ and . This work was performed during the stay of the third author at the Max-Planck-Institute for Mathematics in the Sciences (Leipzig) and also supported by the Russian Fund of Basic Research (grants 05-01-00721, 04-07-90336) and a Priority Research Grant of the Department of Mathematical Sciences of the Russian Academy of Sciences.     

Error-free matrix symmetrizers and equivalent symmetric matrices

A symmetrizer of the matrix A is a symmetric solution X that satisfies the matrix equation XA=AX. An exact matrix symmetrizer is computed by obtaining a general algorithm and superimposing a modified multiple modulus residue arithmetic on this algorithm. A procedure based on computing a symmetrizer to obtain a symmetric matrix, called here an equivalent symmetric matrix, whose eigenvalues are the same as those of a given real nonsymmetric matrix is presented.Supported by CSIR.     

Quadrature rule-based bounds for functions of adjacency matrices

Bounds for entries of matrix functions based on Gauss-type quadrature rules are applied to adjacency matrices associated with graphs. This technique allows to develop inexpensive and accurate upper and lower bounds for certain quantities (Estrada index, subgraph centrality, communicability) that describe properties of networks.     

A direct active set algorithm for large sparse quadratic programs with simple bounds

We show how a direct active set method for solving definite and indefinite quadratic programs with simple bounds can be efficiently implemented for large sparse problems. All of the necessary factorizations can be carried out in a static data structure that is set up before the numeric computation begins. The space required for these factorizations is no larger than that required for a single sparse Cholesky factorization of the Hessian of the quadratic. We propose several improvements to this basic algorithm: a new way to find a search direction in the indefinite case that allows us to free more than one variable at a time and a new heuristic method for finding a starting point. These ideas are motivated by the two-norm trust region problem. Additionally, we also show how projection techniques can be used to add several constraints to the active set at each iteration. Our experimental results show that an algorithm with these improvements runs much faster than the basic algorithm for positive definite problems and finds local minima with lower function values for indefinite problems.Research partially supported by the Applied Mathematical Sciences Research Program (KC-04-02) of the Office of Energy Research of the U.S. Department of Energy under grant DE-FG02-86ER25013.A000.     

The scrambling index of symmetric primitive matrices

A nonnegative square matrix A is primitive if some power A^k>0 (that is, A^k is entrywise positive). The least such k is called the exponent of A. In [2], Akelbek and Kirkland defined the scrambling index of a primitive matrix A, which is the smallest positive integer k such that any two rows of A^k have at least one positive element in a coincident position. In this paper, we give a relation between the scrambling index and the exponent for symmetric primitive matrices, and determine the scrambling index set for the class of symmetric primitive matrices. We also characterize completely the symmetric primitive matrices in this class such that the scrambling index is equal to the maximum value.     

A note on matrices with maximal growth factor for Neville elimination

Neville elimination is a direct method for the solution of linear systems of equations with advantages for some classes of matrices and in the context of pivoting strategies for parallel implementations. The growth factor is an indicator of the numerical stability of an algorithm. In the literature, bounds for the growth factor of Neville elimination with some pivoting strategies have appeared. In this work, we determine all the matrices such that the minimal upper bound of the growth factor of Neville elimination with those pivoting strategies is reached.     

Algorithms for inversion of diagonal plus semiseparable operator matrices

LetH be a Hilbert space andRHH be a bounded linear operator represented by an operator matrix which is a sum of a diagonal and of a semiseparable type of order one operator matrices. We consider three methods for solution of the operator equationRx=y. The obtained results yields new algorithms for solution of integral equations and for inversion of matrices.     

The skew energy of a digraph

We are interested in the energy of the skew-adjacency matrix of a directed graph D, which is simply called the skew energy of D in this paper. Properties of the skew energy of D are studied. In particular, a sharp upper bound for the skew energy of D is derived in terms of the order of D and the maximum degree of its underlying undirected graph. An infinite family of digraphs attaining the maximum skew energy is constructed. Moreover, the skew energy of a directed tree is independent of its orientation, and interestingly it is equal to the energy of the underlying undirected tree. Skew energies of directed cycles under different orientations are also computed. Some open problems are presented.     

Finding minimum height elimination trees for interval graphs in polynomial time

The elimination tree plays an important role in many aspects of sparse matrix factorization. The height of the elimination tree presents a rough, but usually effective, measure of the time needed to perform parallel elimination. Finding orderings that produce low elimination is therefore important. As the problem of finding minimum height elimination tree orderings is NP-hard, it is interesting to concentrate on limited classes of graphs and find minimum height elimination trees for these efficiently. In this paper, we use clique trees to find an efficient algorithm for interval graphs which make an important subclass of chordal graphs. We first illustrate this method through an algorithm that finds minimum height elimination for chordal graphs. This algorithm, although of exponential time complexity, is conceptionally simple and leads to a polynomial-time algorithm for finding minimum height elimination trees for interval graphs.This work was supported by grants from the Norwegian Research Council.