Search

We identify a class of random, dense, $n\times n$ matrices for which LU factorization with any form of pivoting produces a growth factor typically of size at least $n/(4 \log n)$ for large $n$. The...

We propose a two-parameter family of nonsymmetric dense n ⨉ n matrices A(α, β) for which LU factorization without pivoting is numerically stable, and we show how to choose α and β to achieve any value of the ∞-norm condition number. The matrix A(α, β) can be formed from a simple formula in O(n²) flops. The matrix is suitable for use in the HPL-AI Mixed-Precision Benchmark, which requires an extreme scale test matrix (dimension n > 10⁷) that has a controlled condition number and can be safely used in LU factorization without pivoting. It is also of interest as a general-purpose test matrix.

The roundoff-error-free (REF) LU factorization, along with the REF forward and backward substitution algorithms, allows a rational system of linear equations to be solved exactly and efficiently. T...

Motivated by the demand in machine learning, modern computer hardware is increas- ingly supporting reduced precision floating-point arithmetic, which provides advantages in speed, energy, and memory usage over single and double precision. Given the availability of such hardware, mixed precision algorithms that work in single or double precision but carry out part of a compu- tation in half precision are now of great interest for general scientific computing tasks. Because of the limited range of half precision arithmetic, in which positive numbers lie between 6 × 10−8 and 7 × 104, a straightforward rounding of single or double precision data into half precision can lead to overflow, underflow, or subnormal numbers being generated, all of which are undesirable. We develop an algorithm for converting a matrix from single or double precision to half precision. It first applies two-sided diagonal scaling in order to equilibrate the matrix (that is, to ensure that every row and column has ∞-norm 1), then multiplies by a scalar to bring the largest element within a factor θ ≤ 1 of the overflow level, and finally rounds to half precision. The second step ensures that full use is made of the limited range of half precision arithmetic, and θ must be chosen to allow sufficient headroom for subsequent computations. We apply the new algorithm to GMRES-based iterative re- finement (GMRES-IR), which solves a linear system Ax = b with single or double precision data by LU factorizing A in half precision and carrying out iterative refinement with the correction equations solved by GMRES preconditioned with the low precision LU factors. Previous implementations of this algorithm have used a crude conversion to half precision that our experiments show can cause slow convergence of GMRES-IR for badly scaled matrices or failure to converge at all. The new conversion algorithm computes ∞-norms of rows and columns of the matrix and its cost is negligible in the context of LU factorization. We show that it leads to faster convergence of GMRES-IR for badly scaled matrices and thereby allows a much wider class of problems to be solved.

Incomplete factorization is one of the most effective general-purpose preconditioning methods for Krylov subspace solvers for large sparse systems of linear equations. Two techniques for enhancing the robustness and performance of incomplete LU factorization for sparse unsymmetric systems are described. A block incomplete factorization algorithm based on the Crout variation of LU factorization is presented. The algorithm is suitable for incorporating threshold-based dropping as well as unrestricted partial pivoting, and it overcomes several limitations of existing incomplete LU factorization algorithms with and without blocking. It is shown that blocking has a three-pronged impact: it speeds up the computation of incomplete factors and the solution of the associated triangular systems, it permits denser and more robust factors to be computed economically, and it permits a trade-off with the restart parameter of GMRES to further improve the overall speed and robustness. A highly effective heuristic for imp...

Partial integro-differential equation (PIDE) formulations are often preferable for pricing options under models with stochastic volatility and jumps, especially for American-style option contracts. We consider the pricing of options under such models, namely the Bates model and the so-called stochastic volatility with contemporaneous jumps (SVCJ) model. The nonlocality of the jump terms in these models leads to matrices with full matrix blocks. Standard discretization methods are not viable directly since they would require the inversion of such a matrix. Instead, we adopt a two-step implicit-explicit (IMEX) time discretization scheme, the IMEX-CNAB scheme, where the jump term is treated explicitly using the second-order Adams–Bashforth (AB) method, while the rest is treated implicitly using the Crank–Nicolson (CN) method. The resulting linear systems can then be solved directly by employing LU decomposition. Alternatively, the systems can be iterated under a scalable algebraic multigrid (AMG) method. For pricing American options, LU decomposition is employed with an operator splitting method for the early exercise constraint or, alternatively, a projected AMG method can be used to solve the resulting linear complementarity problems. We price European and American options in numerical experiments, which demonstrate the good efficiency of the proposed methods. peerReviewed

We present the LU decomposition with panel rank revealing pivoting (LU_PRRP), an LU factorization algorithm based on strong rank revealing QR panel factorization. LU_PRRP is more stable than Gaussian elimination with partial pivoting (GEPP). Our extensive numerical experiments show that the new factorization scheme is as numerically stable as GEPP in practice, but it is more resistant to pathological cases and easily solves the Wilkinson matrix and the Foster matrix. We also present CALU_PRRP, a communication avoiding version of LU_PRRP that minimizes communication. CALU_PRRP is based on tournament pivoting, with the selection of the pivots at each step of the tournament being performed via strong rank revealing QR factorization. CALU_PRRP is more stable than CALU, the communication avoiding version of GEPP. CALU_PRRP is also more stable in practice and is resistant to pathological cases on which GEPP and CALU fail., Comment: No. RR-7867 (2012)

In 1981 Hong and Kung proved a lower bound on the amount of communication (amount of data moved between a small, fast memory and large, slow memory) needed to perform dense, n-by-n matrix multiplication using the conventional O(n3) algorithm, where the input matrices were too large to fit in the small, fast memory. In 2004 Irony, Toledo, and Tiskin gave a new proof of this result and extended it to the parallel case (where communication means the amount of data moved between processors). In both cases the lower bound may be expressed as Ω(#arithmetic_operations/M), where M is the size of the fast memory (or local memory in the parallel case). Here we generalize these results to a much wider variety of algorithms, including LU factorization, Cholesky factorization, LDLT factorization, QR factorization, the Gram–Schmidt algorithm, and algorithms for eigenvalues and singular values, i.e., essentially all direct methods of linear algebra. The proof works for dense or sparse matrices and for sequential or para...

Quantum Monte Carlo (QMC) methods are often used to calculate properties of many body quantum systems. The main cost of many QMC methods, for example the variational Monte Carlo (VMC) method, is in constructing a sequence of Slater matrices and computing the ratios of determinants for successive Slater matrices. Recent work has improved the scaling of constructing Slater matrices for insulators so that the cost of constructing Slater matrices in these systems is now linear in the number of particles, whereas computing determinant ratios remains cubic in the number of particles. With the long term aim of simulating much larger systems, we improve the scaling of computing the determinant ratios in the VMC method for simulating insulators by using preconditioned iterative solvers. The main contribution of this paper is the development of a method to efficiently compute for the Slater matrices a sequence of preconditioners that make the iterative solver converge rapidly. This involves cheap preconditioner updates, an effective reordering strategy, and a cheap method to monitor instability of ILUTP preconditioners. Using the resulting preconditioned iterative solvers to compute determinant ratios of consecutive Slater matrices reduces the scaling of QMC algorithms from O(n^3) per sweep to roughly O(n^2), where n is the number of particles, and a sweep is a sequence of n steps, each attempting to move a distinct particle. We demonstrate experimentally that we can achieve the improved scaling without increasing statistical errors. Our results show that preconditioned iterative solvers can dramatically reduce the cost of VMC for large(r) systems., 24 pages, 10 figures

We introduce a methodology for obtaining inventories of error results for families of numerical dense linear algebra algorithms. The approach for deriving the analyses is goal-oriented, systematic, and layered. The presentation places the analysis side-by-side with the algorithm so that it is obvious where roundoff error is introduced. The approach supports the analysis of more complex algorithms, such as the blocked LU factorization. For this operation we derive a tighter error bound than has been previously reported in the literature.

This paper presents two methods based on domain decomposition concepts for determining the diagonal of the inverse of specific matrices. The first uses a divide-and-conquer principle and the Sherman-Morrison-Woodbury formula and assumes that the matrix can be decomposed into a $2 \times 2$ block-diagonal matrix and a low-rank matrix. The second method is a standard domain decomposition approach in which local solves are combined with a global correction. Both methods can be successfully combined with iterative solvers and sparse approximation techniques. The efficiency of the methods usually depends on the specific implementation, which should be fine-tuned for different test problems. Preliminary results for some two-dimensional (2D) problems are reported to illustrate the proposed methods.

We discuss and compare two greedy algorithms that compute discrete versions of Fekete-like points for multivariate compact sets by basic tools of numerical linear algebra. The first gives the so-called approximate Fekete points by QR factorization with column pivoting of Vandermonde-like matrices. The second computes discrete Leja points by LU factorization with row pivoting. Moreover, we study the asymptotic distribution of such points when they are extracted from weakly admissible meshes.

In this paper we discuss a hypergraph-based unsymmetric nested dissection (HUND) ordering for reducing the fill-in incurred during Gaussian elimination. It has several important properties. It takes a global perspective of the entire matrix, as opposed to local heuristics. It takes into account the asymmetry of the input matrix by using a hypergraph to represent its structure. It is suitable for performing Gaussian elimination in parallel, with partial pivoting. This is possible because the row permutations performed due to partial pivoting do not destroy the column separators identified by the nested dissection approach. The hypergraph nested dissection approach is essentially equivalent to graph nested dissection on the matrix $A^{T}A$, but we need only the original matrix $A$ and never form the usually denser matrix $A^{T}A$. The usage of hypergraphs in our approach is fairly standard, and HUND can be implemented by calling an existing hypergraph partitioner that uses recursive bisection. Our implementation uses local reordering constrained column approximate minimum degree (CCOLAMD) to further improve the ordering. We also explain how weighted matching (HSL routine MC64) can be used in this context. Experimental results on 27 medium and large size matrices with highly unsymmetric structures compare our approach to four other well-known reordering algorithms. The results show that it provides a robust reordering algorithm, in the sense that it is the best or close to the best (often within 10%) of all the other methods, in particular on matrices with highly unsymmetric structures.

In this paper we consider two structure prediction problems of interest in Gaussian elimination with partial pivoting of sparse matrices. First, we consider the problem of determining the nonzero structure of the factors $L$ and $U$ during the factorization. We present an exact prediction of the structure that identifies some numeric cancellations appearing during Gaussian elimination. The numeric cancellations are related to submatrices of the input matrix $A$ that are structurally singular, that is, singular due to the arrangements of their nonzeros, and independent of their numerical values. Second, we consider the problem of estimating upper bounds for the structure of $L$ and $U$ prior to the numerical factorization. We present tight exact bounds for the nonzero structure of $L$ and $U$ of Gaussian elimination with partial pivoting $PA = LU$ under the assumption that the matrix $A$ satisfies a combinatorial property, namely, the Hall property, and that the nonzero values in $A$ are algebraically independent of each other. This complements existing work showing that a structure called the row merge graph represents a tight bound for the nonzero structure of $L$ and $U$ under a stronger combinatorial assumption, namely, the strong Hall property. We also show that the row merge graph represents a tight symbolic bound for matrices satisfying only the Hall property.

Computing the null space of a sparse matrix is an important part of some computations, such as embeddings and parametrization of meshes. We propose an efficient and reliable method to compute an orthonormal basis of the null space of a sparse square or rectangular matrix (usually with more rows than columns). The main computational component in our method is a sparse $LU$ factorization with partial pivoting of the input matrix; this factorization is significantly cheaper than the $QR$ factorization used in previous methods. The paper analyzes important theoretical aspects of the new method and demonstrates experimentally that it is efficient and reliable.

We describe a fill‐reducing ordering algorithm for sparse, nonsymmetric LU factorizations, where the pivots are restricted to the diagonal and are selected greedily. The ordering algorithm uses only the structural information. Most of the existing methods are based on some type of symmetrization of the original matrix. Our algorithm exploits the nonsymmetric structure of the given matrix as much as possible. The new algorithm is thus more complex than classical symmetric orderings, but we show that our algorithm can be implemented in space bounded by the number of nonzero entries in the original matrix, and has the same time complexity as the analogous algorithms for symmetric matrices. We provide numerical experiments to demonstrate the ordering quality and the runtime of the new ordering algorithm.

We present a family of ordering algorithms that can be used as a preprocessing step prior to performing sparse LU factorization. The ordering algorithms simultaneously achieve the objectives of selecting numerically good pivots and preserving the sparsity. We describe the algorithmic properties and challenges in their implementation. By mixing the two objectives we show that we can reduce the amount of fill-in in the factors and reduce the number of numerical problems during factorization. On a set of large unsymmetric real problems, we obtained the median reductions of 12% in the factorization time, of 13% in the size of the LU factors, of 20% in the number of operations performed during the factorization phase, and of 11% in the memory needed by the multifrontal solver MA41-UNS. A byproduct of this ordering strategy is an incomplete LU-factored matrix that can be used as a preconditioner in an iterative solver.

This paper analyzes dropping strategies in a multilevel incomplete LU decomposition context and presents a few strategies for obtaining related ILUs with enhanced robustness. The analysis shows that the incomplete LU factorization resulting from dropping small entries in Gaussian elimination produces a good preconditioner when the inverses of these factors have norms that are not too large. As a consequence a few strategies are developed whose goal is to achieve this feature. A number of "templates" for enabling implementations of these factorizations are presented. Numerical experiments show that the resulting ILUs offer a good compromise between robustness and efficiency.

The PageRank updating algorithm proposed by Langville and Meyer is a special case of an iterative aggregation/disaggregation (SIAD) method for computing stationary distributions of very large Markov chains. It is designed, in particular, to speed up the determination of PageRank, which is used by the search engine Google in the ranking of web pages. In this paper the convergence, in exact arithmetic, of the SIAD method is analyzed. The SIAD method is expressed as the power method preconditioned by a partial LU factorization. This leads to a simple derivation of the asymptotic convergence rate of the SIAD method. It is known that the power method applied to the Google matrix always converges, and we show that the asymptotic convergence rate of the SIAD method is at least as good as that of the power method. Furthermore, by exploiting the hyperlink structure of the web it can be shown that the asymptotic convergence rate of the SIAD method applied to the Google matrix can be made strictly faster than that of the power method.

We revise the dual projective pivot algorithm using sparse rectangular LU factors. In each iteration, the proposed algorithm solves at most three triangular systems, while simplex algorithms solve four. Instead of the classical basis, it uses a so-called pseudobasis (a rectangular matrix having fewer columns than rows), thereby solving smaller linear systems with a potentially improved stability compared to simplex algorithms. Most importantly, it generates good search directions at a low cost. We report encouraging computational results on a set of 50 Netlib standard test problems as well as a set of 15 much larger real-world problems. A code named RDPPA 1.10 based on the proposed algorithm outperformed MINOS 5.51 significantly in terms of both iterations and run time. In particular, it appears that a high proportion of degenerate iterations need not imply many total iterations (contradicting the common belief).

Component- and normwise perturbation bounds for the block LU factorization and block LDL$^*$ factorization of Hermitian matrices are presented. We also obtain, as a consequence, perturbation bounds for the usual pointwise LU, LDL$^*$, and Cholesky factorizations. Some of these latter bounds are already known, but others improve previous results. All the bounds presented are easily proved by using series expansions. Given a square matrix $A=L U$ having the LU factorization, and a perturbation $E$, the LU factors of the matrix $A+E= \widetilde{L} \widetilde{U}$ are written as two convergent series of matrices: $\widetilde{L} = \sum_{k=0}^{\infty} L_k$ and $\widetilde{U} = \sum_{k=0}^{\infty} U_k$, where $L_k = O(\|E\|^k)$, $U_k = O(\|E\|^k)$, and $L_0 = L$, $U_0 = U$. We present expressions for the matrices $L_k$ and $U_k$ in terms of $L$, $U$, and $E$. The domain and the rate of convergence of these series are studied. Simple bounds on the remainders of any order of these series are found, which significantly improve the bounds on the second-order terms existing in the literature. This is useful when first-order perturbation analysis is used.

This paper presents a preconditioning method based on combining two-sided permutations with a multilevel approach. The nonsymmetric permutation exploits a greedy strategy to put large entries of the matrix in the diagonal of the upper leading submatrix. The method can be regarded as a complete pivoting version of the incomplete LU factorization. This leads to an effective incomplete factorization preconditioner for general nonsymmetric, irregularly structured, sparse linear systems.

In this paper we present a new family of methods for the effective numerical solution of the incompressible unsteady Navier--Stokes equations. These methods resort to an algebraic splitting of the discretized problem based on inexact LU block factorizations of the corresponding matrix (following [A. Quarteroni, F. Saleri, and A. Veneziani, Comput. Methods Appl. Mech. Engrg., 188 (2000), pp. 505--526]. In particular, we will start from inexact algebraic factorizations of algebraic Chorin--Temam and Yosida type and introduce a pressure correction step aimed at improving the time accuracy. One of the schemes obtained in this way (the algebraic Chorin--Temam pressure correction method) resembles a method previously introduced in the framework of differential projection schemes (see [L. Timmermans, P. Minev, and F. V. de Vosse, Internat. J. Numer. Methods Fluids, 22 (1996), pp. 673--688], [A. Prohl, Projection and Quasi-Compressibility Methods for Solving the Incompressible Navier--Stokes Equations, Teubner, Stuttgart, 1997]. The stability and the dependence of splitting error on the time step of the new methods is investigated and tested on several numerical cases.

The elimination tree of a symmetric matrix plays an important role in sparse matrix factorization. By using paths instead of edges to define the tree, we generalize this structure to unsymmetric matrices while retaining many of its properties. If we use a tree traversal to reorder a matrix into a bordered block triangular form, the structure has further desirable properties relevant to a sparse LU factorization of the reordered matrix. When pivoting is required for stability, the tree changes only locally if the choice of pivot is suitably restricted.

In this paper we present a partially orthogonal decomposition for a matrix A. Using this decomposition the linear least squares problem is reduced to solving two linear systems. The matrix of the first system is symmetric and positive definite, and the matrix of the second system is nonsingular upper triangular. We show that this approach can provide computational savings.

We investigate the problem of permuting a sparse rectangular matrix into block-diagonal form. Block-diagonal form of a matrix grants an inherent parallelism for solving the deriving problem, as recently investigated in the context of mathematical programming, LU factorization, and QR factorization. To represent the nonzero structure of a matrix, we propose bipartite graph and hypergraph models that reduce the permutation problem to those of graph partitioning by vertex separator and hypergraph partitioning, respectively. Our experiments on a wide range of matrices, using the state-of-the-art graph and hypergraph partitioning tools MeTiS and PaToH\@, revealed that the proposed methods yield very effective solutions both in terms of solution quality and runtime.

We present a new method, called reversed double ordering, for reordering arbitrary matrices prior to LU factorization. This reordering creates a variable-band matrix. We compare the fill-in of the LU factorization for sparse matrices with respect to reversed double ordering, column minimum degree ordering, and the reversed Cuthill--McKee algorithm. Moreover, we combine the first two reorderings with good success.

In this paper, a general-purpose block LU preconditioner for saddle point problems is presented. The main difference between the approach presented here and that of other studies is that an explicit, accurate approximation of the Schur complement matrix is efficiently computed. This is used to compute a preconditioner to the Schur complement matrix, which in turn defines a preconditioner for the global system. A number of different variants are developed and results are reported for a few linear systems arising from CFD applications.

This short communication considers the LU factorization with partial pivoting and shows that an all-at-once result is possible for the structure prediction of the column dependencies in L and U. Specifically, we prove that for every square strong Hall matrix A there exists a permutation P such that every edge of its column elimination tree corresponds to a symbolic nonzero in the upper triangular factor U. In the symbolic sense, this resolves a conjecture of Gilbert and Ng [Graph Theory and Sparse Matrix Computation, A. George, J. R. Gilbert, and J. W. H. Liu, eds., Springer-Verlag, New York, 1994, pp. 107--139].

Sparse matrices which arise in many applications often possess a block structure that can be exploited in iterative and direct solution methods. These block-matrices have as their entries small dense blocks with constant or variable dimensions. Block versions of incomplete LU factorizations which have been developed to take advantage of such structures give rise to a class of preconditioners that are among the most effective available. This paper presents general techniques for automatically determining block structures in sparse matrices. A standard "graph compression" algorithm used in direct sparse matrix methods is considered along with two other algorithms which are also capable of unraveling approximate block structures.

A new solution technique is proposed for linear systems with large dense matrices of a certain class including those that come from typical integral equations of potential theory. This technique combines Kronecker product approximation and wavelet sparsification for the Kronecker product factors. The user is only required to supply a procedure for computation of each entry of the given matrix. The main sources of efficiency are the incomplete cross approximation procedure adapted from the mosaic-skeleton method of the second author and data-sparse preconditioners (the incomplete LU decomposition with dynamic choice of the fill-in structure with a prescribed threshold and the inverse Kronecker product preconditioner) constructed for the sum of Kronecker products of sparsified finger-like matrices computed by the discrete wavelet transform. In some model, but quite representative, examples the new technique allowed us to solve dense systems with more than 1 million unknowns in a few minutes on a personal computer with 1 Gbyte operative memory.

This paper discusses some relationships between ILU factorization techniques and factored sparse approximate inverse techniques. While ILU factorizations compute approximate LU factors of the coefficient matrix A, approximate inverse techniques aim at building triangular matrices Z and W such that $W^\top AZ$ is approximately diagonal. The paper shows that certain forms of approximate inverse techniques amount to approximately inverting the triangular factors obtained from some variants of ILU factorization of the original matrix. A few useful applications of these relationships will be discussed.

We describe the implementation and performance of a novel fill-minimization ordering technique for sparse LU factorization with partial pivoting. The technique was proposed by Gilbert and Schreiber in 1980 but never implemented and tested. Like other techniques for ordering sparse matrices for LU with partial pivoting, our new method preorders the columns of the matrix (the row permutation is chosen by the pivoting sequence during the numerical factorization). Also like other methods, the column permutation Q that we select is a permutation that attempts to reduce the fill in the Cholesky factor of QT ATAQ. Unlike existing column-ordering techniques, which all rely on minimum-degree heuristics, our new method is based on a nested-dissection ordering of A T A. Our algorithm, however, never computes a representation of AT A, which can be expensive. We only work with a representation of A itself. Our experiments demonstrate that the method is efficient and that it can reduce fill significantly relative to the best existing methods. The method reduces the LU running time on some very large matrices (tens of millions of nonzeros in the factors) by more than a factor of 2.

A well-known approach to compute the LU factorization of a general unsymmetric matrix A is to build the elimination tree associated with the pattern of the symmetric matrix A + A rm T and use it as a computational graph to drive the numerical factorization. This approach, although very efficient on a large range of unsymmetric matrices, does not capture the unsymmetric structure of the matrices. We introduce a new algorithm which detects and exploits the structural unsymmetry of the submatrices involved during the process of the elimination tree. We show that with the new algorithm significant gains both in memory and in time to perform the factorization can be obtained.

We describe an algebraic multilevel multigraph algorithm. Many of the multilevel components are generalizations of algorithms originally applied to general sparse Gaussian elimination. Indeed, general sparse Gaussian elimination with minimum degree ordering is a limiting case of our algorithm. Our goal is to develop a procedure which has the robustness and simplicity of use of sparse direct methods, yet offers the opportunity to obtain the optimal or near-optimal complexity typical of classical multigrid methods.

In this paper, we analyze the main features and discuss the tuning of the algorithms for the direct solution of sparse linear systems on distributed memory computers developed in the context of a long term European research project. The algorithms use a multifrontal approach and are especially designed to cover a large class of problems. The problems can be symmetric positive definite, general symmetric, or unsymmetric matrices, both possibly rank deficient, and they can be provided by the user in several formats. The algorithms achieve high performance by exploiting parallelism coming from the sparsity in the problem and that available for dense matrices. The algorithms use a dynamic distributed task scheduling technique to accommodate numerical pivoting and to allow the migration of computational tasks to lightly loaded processors. Large computational tasks are divided into subtasks to enhance parallelism. Asynchronous communication is used throughout the solution process to efficiently overlap communication with computation. We illustrate our design choices by experimental results obtained on an SGI Origin 2000 and an IBM SP2 for test matrices provided by industrial partners in the PARASOL project.

We describe a parallel algorithm for computing incomplete factor (ILU) preconditioners. The algorithm attains a high degree of parallelism through graph partitioning and a two-level ordering strategy. Both the subdomains and the nodes within each subdomain are ordered to preserve concurrency. We show through an algorithmic analysis and through computational results that this algorithm is scalable. Experimental results include timings on three parallel platforms for problems with up to 20 million unknowns running on up to 216 processors. The resulting preconditioned Krylov solvers have the desirable property that the number of iterations required for convergence is insensitive to the number of processors.

We present a roundoff error analysis of a null space method for solving quadratic programming minimization problems. This method combines the use of a direct LU factorization of the constraints with an iterative solver on the corresponding null space. Numerical experiments are presented which give evidence of the good performance of the algorithm on sparse matrices.

We present a family of algorithms for computing symmetric rank-revealing VSV decompositions based on triangular factorization of the matrix. The VSV decomposition consists of a middle symmetric matrix that reveals the numerical rank in having three blocks with small norm, plus an orthogonal matrix whose columns span approximations to the numerical range and null space. We show that for semidefinite matrices the VSV decomposition should be computed via the ULV decomposition, while for indefinite matrices it must be computed via a URV-like decomposition that involves hypernormal rotations.

A refinement of the polar decomposition of a nonsingular matrix A is considered. Here A is written as a product of unitary U, Hermitian and positive definite P which has unit diagonal, and diagonal positive D. It is shown that such a decomposition exists and is unique. Rectangular and singular cases are also considered. Then a simple fixed point iteration using SVD is given to compute this decomposition. Also, implementation of the Newton's method is discussed. The refined polar decomposition can be used to parameterize the orbit of a matrix with distinct eigenvalues.

We introduce a new preconditioner for elliptic PDEs on unstructured meshes. Using a wavelet-inspired basis we compress the inverse of the matrix, allowing an effective sparse approximate inverse by solving the sparsity vs. accuracy conflict. The key issue in this compression is to use second generation wavelets which can be adapted to the unstructured mesh, the true boundary conditions, and even the PDE coefficients. We also show how this gives a new perspective on multiresolution algorithms such as multigrid, interpreting the new preconditioner as a variation on node-nested multigrid. In particular, we hope the new preconditioner will combine the best of both worlds: fast convergence when multilevel methods can succeed but with robust performance for more difficult problems. The rest of the paper discusses the core issues for the preconditioner: ordering and construction of a factored approximate inverse in the multiresolution basis, robust interpolation on unstructured meshes, automatic mesh coarsening, and purely algebraic alternatives. Some exploratory numerical experiments suggest the superiority of the new basis over the standard basis for several tough problems, including discontinuous anisotropic coefficients, strong convection, and indefinite reaction problems on unstructured meshes, with scalability like hierarchical basis methods achieved.

Static symbolic factorization coupled with supernode partitioning and asynchronous computation scheduling can achieve high gigaflop rates for parallel sparse LU factorization with partial pivoting. This paper studies properties of elimination forests and uses them to optimize supernode partitioning/amalgamation and execution scheduling. It also proposes supernodal matrix multiplication to speed up kernel computation by retaining the BLAS-3 level efficiency and avoiding unnecessary arithmetic operations. The experiments show that our new design with proper space optimization, called S+, improves our previous solution substantially and can achieve up to 10 GFLOPS on 128 Cray T3E 450MHz nodes.

This paper presents a new algorithm for accurate floating-point computation of the singular value decomposition (SVD) of the product $A = B^{\tau} S C,$ where $ B\in {\hbox{{\bf R}$^{p\times m} $}},$ $C\in {\hbox{{\bf R}$^{q\times n}$}},$ $S\in {\hbox{{\bf R}$^{p\times q}$}},$ and $p\leq m,$ $q\leq n$.\ The new algorithm uses diagonal scalings, the LU factorization with complete pivoting, the QR factorization with column pivoting, and matrix multiplication to replace $A$ by $A' = B'^{\tau}S'C',$ where A and A' have the same singular values and the matrix A' is computed explicitly. The singular values of A' are computed using the Jacobi SVD algorithm. It is shown that the accuracy of the new algorithm is determined by (i) the accuracy of the QR factorizations of $B^{\tau}$ and $C^{\tau}$; (ii) the accuracy of the LU factorization with complete pivoting of S; and (iii) the accuracy of the computation of the SVD of a matrix $A'$ with moderate $\min_{D=\diag} \kappa_2(A'D)$.\ Theoretical analysis and numerical evidence show that, in the case of rank (B)= rank(C)=p and full rank S, the accuracy of the new algorithm is unaffected by replacing B, S, C with, respectively, D1 B, D2SD3,D4 C, where Di, i=1, . . .,4, are arbitrary diagonal matrices. As an application, the paper proposes new accurate algorithms for computing the (H,K)--SVD and (H-1,K)--SVD of S.

We discuss a multilinear generalization of the singular value decomposition. There is a strong analogy between several properties of the matrix and the higher-order tensor decomposition; uniqueness, link with the matrix eigenvalue decomposition, first-order perturbation effects, etc., are analyzed. We investigate how tensor symmetries a ect the decomposition and propose a multilinear generalization of the symmetric eigenvalue decomposition for pair-wise symmetric tensors. ispartof: SIAM Journal on Matrix Analysis and Applications vol:21 issue:4 pages:1253-1278 status: published

We introduce block versions of the multielimination incomplete LU (ILUM) factorization preconditioning technique for solving general sparse unstructured linear systems. These preconditioners have a multilevel structure and, for certain types of problems, may exhibit properties that are typically enjoyed by multigrid methods. Several heuristic strategies for forming blocks of independent sets are introduced and their relative merits are discussed. The advantages of block ILUM over point ILUM include increased robustness and efficiency. We compare several versions of the block ILUM, point ILUM, and the dual-threshold-based ILUT preconditioners. In particular, tests with some convection-diffusion problems show that it may be possible to obtain convergence that is nearly independent of the Reynolds number as well as of the grid size.

We present a new family of multigraph algorithms, ILU-MG, based upon an incomplete sparse matrix factorization using a particular ordering and allowing a limited amount of fill-in. While much of the motivation for multigraph comes from multigrid ideas, ILU-MG is distinctly different from algebraic multilevel methods. The graph of the sparse matrix A is recursively coarsened by eliminating vertices using a graph model similar to Gaussian elimination. Incomplete factorizations are obtained by allowing only the fill-in generated by the vertex parents associated with each vertex. Multigraph is numerically compared with algebraic multigrid on some examples arising from discretizations of partial differential equations on unstructured grids.

We investigate several ways to improve the performance of sparse LU factorization with partial pivoting, as used to solve unsymmetric linear systems. We introduce the notion of unsymmetric supernodes to perform most of the numerical computation in dense matrix kernels. We introduce unsymmetric supernode-panel updates and two-dimensional data partitioning to better exploit the memory hierarchy. We use Gilbert and Peierls's depth-first search with Eisenstat and Liu's symmetric structural reductions to speed up symbolic factorization. We have developed a sparse LU code using all these ideas. We present experiments demonstrating that it is significantly faster than earlier partial pivoting codes. We also compare its performance with UMFPACK, which uses a multifrontal approach; our code is very competitive in time and storage requirements, especially for large problems.

In this paper we introduce a new decomposition called the pivoted QLP decomposition. It is computed by applying pivoted orthogonal triangularization to the columns of the matrix X in question to get an upper triangular factor R and then applying the same procedure to the rows of R to get a lower triangular matrix L. The diagonal elements of R are called the R-values of X; those of L are called the L-values. Numerical examples show that the L-values track the singular values of X with considerable fidelity---far better than the R-values. At a gap in the L-values the decomposition provides orthonormal bases of analogues of row, column, and null spaces provided of X. The decomposition requires no more than twice the work required for a pivoted QR decomposition. The computation of R and L can be interleaved, so that the computation can be terminated at any suitable point, which makes the decomposition especially suitable for low-rank determination problems. The interleaved algorithm also suggests a new, efficient 2-norm estimator.

It is shown that maximal growth for Gaussian elimination with partial pivoting as measured in the 2-norm is achieved by orthogonal matrices. A precise bound on that growth is given.

For a real nonsingular n-by-n matrix A, there exists a decomposition $A=V\Pi U$, where $\Pi$ is a permutation matrix and V,U are upper triangular matrices. When $\Pi^TV\Pi$ is lower triangular and U is normalized, such a decomposition is called the left Bruhat decomposition of A. An algorithm for computing the left Bruhat decomposition is given. For classes of matrices introduced by Wilkinson and recently (from a practical application) by Foster that have an exponential growth factor when Gaussian elimination with partial pivoting (GEPP) is applied, left Bruhat decomposition has at most linear growth. A partial pivoting strategy for Bruhat decomposition is also developed, and an explicit equivalence between GEPP and Bruhat decomposition with partial pivoting (BDPP) is derived. This equivalence implies that the growth factor for GEPP on A equals the growth factor for BDPP on $\rho A^T$, where $\rho$ is the permutation matrix that reverses the rows of $A^T$. BDPP is shown to give a growth factor of at most 2 when applied to any matrix for which GEPP gives the maximal growth factor of 2 n - 1.