Pivoting Research Articles

Undulant-block elimination and integer-preserving matrix inversion

A new formulation for LU decomposition allows efficient representation of intermediate matrices while eliminating blocks of various sizes, i.e. during “undulant-block” elimination. Its efficiency arises from its design for block encapsulization, implicit in data structures that are convenient both for process scheduling and for memory management. Row/column permutations that can destroy such encapsulizations are deferred. Its algorithms, expressed naturally as functional programs, are well suited to parallel and distributed processing. A given matrix A is decomposed into two matrices (in the space of just one), plus two permutations. The permutations, P and Q, are the row/column rearrangements usual to complete pivoting. The principal results are L and U′, where L is properly lower quasi-triangular; U′ is upper quasi-triangular with its quasi-diagonal being the inverse of that of U from the usual factorization ( PAQ = ( I − L) U), and its proper upper portion identical to U. The matrix result is L + U′. Algorithms for solving linear systems and matrix inversion follow directly. An example of a motivating data structure, the quadtree representation for matrices, is reviewed. Candidate pivots for Gaussian elimination under that structure are the subtrees, both constraining and assisting the pivot search, as well as decomposing to independent block/tree operations. The elementary algorithms are provided, coded in Haskell. Finally, an integer-preserving version is presented replacing Bareiss's algorithm with a parallel equivalent. The decomposition of an integer matrix A to integer matrices L ̄ , U ̄ ′ , and d = det A follows L + U′ decomposition, but the follow-on algorithm to compute dA −1 is complicated by the requirement to maintain minimal denominators at every step and to avoid divisions, restricting them to necessarily exact ones.

On the sensitivity of the LU factorization

This paper gives sensitivity analyses by two approaches forL andU in the factorizationA=LU for general perturbations inA which are sufficiently small in norm. By the matrix-vector equation approach, we derive the condition numbers for theL andU factors. By the matrix equation approach we derive corresponding condition estimates. We show how partial pivoting and complete pivoting affect the sensitivity of the LU factorization.

Three-Dimensional Solution for the Diffraction Problem of a Submerged Body in Waves

The paper presents a panel method to evaluate the diffraction problem on a submerged body moving beneath the free surface in regular waves. In the process of numerical calculation, the 1/r term in the potential function is treated by the method submitted by Hess & Smith (1964) with three ways for the integration, and the pulsating source terms due to the free surface effect are solved by using four different series expansions and an integral representation submitted by Telste & Noblesse (1986). With such derivation, the numerical computing techniques become available to calculate the wave exciting force for a submerged body advancing in regular waves at constant speed. Besides, the effect of the steady potential is also considered in the paper. The replacement of the double integration by a single integral in terms of a complex exponential integral is made in the derivative of the steady potential. The source strength on panels is solved efficiently by applying the Gaussian elimination method incorporated with the partial pivoting technique. Results obtained by the present technique are compared with the other methods and good agreement is found. The effect of steady flow is generally small, but it still cannot be neglected in some cases, especially for the slender body in oblique waves. It is also found that using the present technique incorporated with a suitable body mesh distribution can lead to more satisfactory results.

Large‐scale DAE optimization using a simultaneous NLP formulation

AbstractThe differential‐algebraic equation (DAE) optimization problem is transformed to a nonlinear programming problem by applying collocation on finite elements. The resulting problem is solved using a reduced space successive quadratic programming (rSQP) algorithm. Here, the variable space is partitioned into range and null spaces. Partitioning by choosing a pivot sequence for an LU factorization with partial pivoting allows us to detect unstable modes in the DAE system, which can now be stabilized without imposing new boundary conditions. As a result, the range space is decomposed in a single step by exploiting the overall sparsity of the collocation matrix; which performs better than the two‐step condensation method used in standard collocation solvers. To deal with ill‐conditioned constraints, we also extend the rSQP algorithm to include dogleg steps for the range space step that solves the collocation equations. The performance of this algorithm was tested on two well known unstable problems and on three chemical engineering examples including two reactive distillation columns and a plug‐flow reactor with free radicals. One of these is u batch column where an equilibrium reaction takes place. The second reactive distillation problem is the startup of a continuous column with competitive reactions. These optimization problems, which include more than 150 DAEs, ure solved in less than 7 CPU minutes on workstation class computers.

Stable and Efficient Algorithms for Structured Systems of Linear Equations

Recent research shows that structured matrices such as Toeplitz and Hankel matrices can be transformed into a different class of structured matrices called Cauchy-like matrices using the FFTor other trigonometric transforms. Gohberg, Kailath, and Olshevsky [Math. Comp., 64 (1995), pp. 1557--1576] demonstrate numerically that their fast variation of the straightforward Gaussian elimination with partial pivoting (GEPP) procedure on Cauchy-like matrices is numerically stable. Sweet and Brent [Adv. Signal Proc. Algorithms, 2363 (1995), pp. 266--280] show that the error growth in this variation could be much larger than would be encountered with straightforward GEPP in certain cases. In this paper, we present a modified algorithm that avoids such extra error growth and can perform a fast variation of Gaussian elimination with complete pivoting (GECP). Our analysis shows that it is both efficient and numerically stable, provided that the element growth in the computed factorization is not large. We also present a more efficient variation of this algorithm and discuss implementation techniques that further reduce execution time. Our numerical experiments show that this variation is highly efficient and numerically stable.

Efficient sparse LU factorization with partial pivoting on distributed memory architectures

A sparse LU factorization based on Gaussian elimination with partial pivoting (GEPP) is important to many scientific applications, but it is still an open problem to develop a high performance GEPP code on distributed memory machines. The main difficulty is that partial pivoting operations dynamically change computation and nonzero fill-in structures during the elimination process. This paper presents an approach called S* for parallelizing this problem on distributed memory machines. The S* approach adopts static symbolic factorization to avoid run-time control overhead, incorporates 2D L/U supemode partitioning and amalgamation strategies to improve caching performance, and exploits irregular task parallelism embedded in sparse LU using asynchronous computation scheduling. The paper discusses and compares the algorithms using 1D and 2D data mapping schemes, and presents experimental studies on Cray-T3D and T3E. The performance results for a set of nonsymmetric benchmark matrices are very encouraging, and S* has achieved up to 6.878 GFLOPS on 128 T3E nodes. To the best of our knowledge, this is the highest performance ever achieved for this challenging problem and the previous record was 2.583 GFLOPS on shared memory machines.

Bruhat Decomposition and Numerical Stability

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.

The growth factor and efficiency of Gaussian elimination with rook pivoting

Gaussian elimination is among the most widely used tools in scientific computing. Gaussian elimination with partial pivoting requires only O( n 2) comparisons beyond the work required in Gaussian elimination with no pivoting but can, in principle, have error growth that is exponential in the matrix size n. Gaussian elimination with complete pivoting, on the other hand, cannot have exponential error growth but requires O( n 3) comparisons beyond the work required by Gaussian elimination with no pivoting. Numerical experiments suggest that Gaussian elimination with rook pivoting is between partial pivoting and complete pivoting in terms of efficiency and accuracy. In this paper we prove that rook pivoting cannot have exponential error growth. We also introduce a combination of partial pivoting and rock pivoting which we call Gaussian elimination with partial rook pivoting and we prove the partial rook pivoting cannot have exponential error growth. We include numerical experiments showing that on a serial computer the run times for cook pivoting are almost always close to those of partial pivoting and the run times for partial rook pivoting appear to be the same as those of partial pivoting.

Recursion leads to automatic variable blocking for dense linear-algebra algorithms

We describe some modifications of the LAPACK dense linear-algebra algorithms using recursion. Recursion leads to automatic variable blocking. LAPACK's level-2 versions transform into level-3 codes by using recursion. The new recursive codes are written in FORTRAN 77, which does not support recursion as a language feature. Gaussian elimination with partial pivoting and Cholesky factorization are considered. Very clear algorithms emerge with the use of recursion. The recursive codes do exactly the same computation as the LAPACK codes, and a single recursive code replaces both the level-2 and level-3 versions of the corresponding LAPACK codes. We present an analysis of the recursive algorithm in terms of both FLOP count and storage usage. The matrix operands are more “squarish” using recursion. The total area of the submatrices used in the recursive algorithm is less than the total area used by the LAPACK level-3 right-/left-looking algorithms. We quantify the difference; we also quantify how the FLOPS are computed. Also, we show that the algorithms exhibit high-performance on RISC-type processors. In fact, except for small matrices, the recursive version outperforms the level-3 LAPACK versions of DGETRF and DPOTRF on an RS/6000™ workstation. For the level-2 versions, the performance gain approaches a factor of 3. We also demonstrate that a change to the LAPACK DLASWP routine can improve the performance of both the recursive version and DGETRF by more than 15 percent.

Locality of Reference in LU Decomposition with Partial Pivoting

This paper presents a new partitioned algorithm for LU decomposition with partial pivoting. The new algorithm, called the recursively partitioned algorithm, is based on a recursive partitioning of the matrix. The paper analyzes the locality of reference in the new algorithm and the locality of reference in a known and widely used partitioned algorithm for LU decomposition called the right-looking algorithm. The analysis reveals that the new algorithm performs a factor of $\Theta(\sqrt{M/n})$ fewer I/O operations (or cache misses) than the right-looking algorithm, where $n$ is the order of the matrix and $M$ is the size of primary memory. The analysis also determines the optimal block size for the right-looking algorithm. Experimental comparisons between the new algorithm and the right-looking algorithm show that an implementation of the new algorithm outperforms a similarly coded right-looking algorithm on six different RISC architectures, that the new algorithm performs fewer cache misses than any other algorithm tested, and that it benefits more from Strassen's matrix-multiplication algorithm.

Stability of the Gauss-Huard algorithm with partial pivoting

This paper considers elimination methods to solve dense linear systems, in particular a variant of Gaussian elimination due to Huard [13]. This variant reduces the system to an equivalent diagonal system just like Gauss-Jordan elimination, but does not require more floating-point operations than Gaussian elimination. To preserve stability, a pivoting strategy using column interchanges, proposed by Hoffmann [10], is incorporated in the original algorithm. An error analysis is given showing that Huard’s elimination method is as stable as Gauss-Jordan elimination with the appropriate pivoting strategy. This result is proven in a similar way as the proof of stability for Gauss-Jordan elimination given in [4]. Numerical experiments are reported which verify the theoretical error analysis of the Gauss-Huard algorithm.

Space and time efficient execution of parallel irregular computations

Solving problems of large sizes is an important goal for parallel machines with multiple CPU and memory resources. In this paper, issues of efficient execution of overhead-sensitive parallel irregular computation under memory constraints are addressed. The irregular parallelism is modeled by task dependence graphs with mixed granularities. The trade-off in achieving both time and space efficiency is investigated. The main difficulty of designing efficient run-time system support is caused by the use of fast communication primitives available on modern parallel architectures. A run-time active memory management scheme and new scheduling techniques are proposed to improve memory utilization while retaining good time efficiency, and a theoretical analysis on correctness and performance is provided. This work is implemented in the context of RAPID system [5] which provides run-time support for parallelizing irregular code on distributed memory machines and the effectiveness of the proposed techniques is verified on sparse Cholesky and LU factorization with partial pivoting. The experimental results on Cray-T3D show that solvable problem sizes can be increased substantially under limited memory capacities and the loss of execution efficiency caused by the extra memory managing overhead is reasonable.

Probabilistic Analysis of Gaussian Elimination Without Pivoting

The numerical instability of Gaussian elimination is proportional to the size of the $L$ and $U$ factors that it produces. The worst-case bounds are well known. For the case without pivoting, breakdowns can occur and it is not possible to provide {a priori} bounds for $L$ and $U$. For the partial pivoting case, the worst-case bound is $O(2^m)$, where $m$ is the size of the system. Yet these worst-case bounds are seldom achieved, and in particular Gaussian elimination with partial pivoting is extremely stable in practice. Surprisingly, there has been relatively little theoretical study of the ``average'' case behavior. The purpose of our paper is to provide a probabilistic analysis of the case without pivoting. The distribution we use for the entries of $A$ is the normal distribution with mean $0$ and unit variance. We first derive the distributions of the entries of $L$ and $U$. Based on this, we prove that the probability of the occurrence of a pivot less than $\epsilon$ in magnitude is $O(\epsilon)$. We also prove that the probabilities ${\rm Prob}(||U||_\infty/||A||_\infty > m^{2.5})$ and ${\rm Prob}(||L||_\infty > m^{3})$ decay algebraically to zero as $m$ tends to infinity. Numerical experiments are presented to support the theoretical results.

Diagonal pivoting for partially reconstructible Cauchy-like matrices, with applications to Toeplitz-like linear equations and to boundary rational matrix interpolation problems

In an earlier paper we exploited the displacement structure of Cauchy-like matrices to derive for them a fast O( n 2) implementation of Gaussian elimination with partial pivoting. One application is to the rapid and numerically accurate solution of linear systems with Toeplitz-like coefficient matrices, based on the fact that the latter can be transformed into Cauchy-like matrices by using the fast Fourier, sine, or cosine transform. However, symmetry is lost in the process, and the algorithm given is not optimal for Hermitian coefficient matrices. In this paper we present a new fast O( n 2) implementation of symmetric Gaussian elimination with partial diagonal pivoting for Hermitian Cauchy-like matrices, and show how to transform Hermitian Toeplitz-like matrices to Hermitian Cauchy-like matrices, obtaining algorithms that are twice as fast as those in the earlier work. Numerical experiments indicate that in order to obtain not only fast but also numerically accurate methods, it is advantageous to explore the important case in which the corresponding displacement operators have nontrivial kernels; this situation gives rise to what we call partially reconstructible matrices, which are introduced and studied in the present paper. We extend the transformation technique and the generalized Schur algorithms (i.e., fast displacement-based implementations of Gaussian elimination) to partially reconstructible matrices. We show by a variety of computed examples that the incorporation of diagonal pivoting methods leads to high accuracy. We focused in this paper on the design of new numerically reliable algorithms for Hermitian Toeplitz-like matrices. However, the proposed algorithms have other important applications; in particular, we briefly describe how they recursively solve a boundary interpolation problem for J-unitary rational matrix functions.

Cholesky factorization of semidefinite Toeplitz matrices

It can be shown directly from consideration of the Schur algorithm that any n × n semidefinite rank r Toeplitz matrix, T, has a factorization T = C r C T r with C r = C 11 C 12 0 0 where C 11 is r × r and upper triangular. This paper explores the reliability of computing such a decomposition with O( nr) complexity using the Schur algorithm and truncating the Cholesky factor after computing the first r rows. The theoretical conclusion is that, as with Gaussian elimination with partial pivoting, the algorithm is stable in a qualified sense: the backward error bounds can be expressed purely in terms of n and r and are independent of the condition number of the leading r × r submatrix T 11 of T, but they grow exponentially in r. Thus, the algorithm is completely reliable when r is small, but assessing its reliability for larger r requires making a judgment as to how realistic the error bounds are in practice. Experimental results are presented which suggest that the bounds are pessimistic, and the algorithm is stable for many types of semidefinite Toeplitz matrices. The analysis and conclusions are similar to those in a paper by Higham dealing with general semidefinite matrices.

Data and Workload Distribution in a Multithreaded Architecture

Matching data distribution to workload distribution is important in improving the performance of distributed-memory multiprocessors. While data and workload distribution can be tailored to fit a particular problem to a particular distributed-memory architecture, it is often difficult to do so for various reasons including complexity of address computation, runtime data movement, and irregular resource usage. This report presents our study on multithreading for distributed-memory multiprocessors. Specifically, we investigate the effects of multithreading ondatadistribution andworkloaddistribution withvariable, thread granularity. Various types of workload distribution strategies are defined along with thread granularity. Several types of data distribution strategies are investigated. These include row-wise cyclic,k-way partial-row cyclic, and blocked distribution. To investigate the performance of multithreading, two problems are selected: highly sequential Gaussian elimination with partial pivoting and highly parallel matrix multiplication. Execution results on the 80-processor EM-4 distributed-memory multiprocessor indicate that multithreading can off set the loss due to the mismatch between data distribution and workload distribution even for sequential and irregular problems while giving high absolute performance.

Stability of the Diagonal Pivoting Method with Partial Pivoting

LAPACK and LINPACK both solve symmetric indefinite linear systems using the diagonal pivoting method with the partial pivoting strategy of Bunch and Kaufman [Math. Comp., 31 (1977), pp. 163--179]. No proof of the stability of this method has appeared in the literature. It is tempting to argue that the diagonal pivoting method is stable for a given pivoting strategy if the growth factor is small. We show that this argument is false in general and give a sufficient condition for stability. This condition is not satisfied by the partial pivoting strategy because the multipliers are unbounded. Nevertheless, using a more specific approach we are able to prove the stability of partial pivoting, thereby filling a gap in the body of theory supporting LAPACK and LINPACK.

Multi-Terminal Binary Decision Diagrams: An Efficient DataStructure for Matrix Representation

In this paper, we discuss the use of binary decision diagrams to represent general matrices. We demonstrate that binary decision diagrams are an efficient representation for every special-case matrix in common use, notably sparse matrices. In particular, we demonstrate that for any matrix, the BDD representation can be no larger than the corresponding sparse-matrix representation. Further, the BDD representation is often smaller than any other conventional special-case representation: for the n×n Walsh matrix, for example, the BDD representation is of size O(log n). No other special-case representation in common use represents this matrix in space less than O(n²). We describe termwise, row, column, block, and diagonal selection over these matrices, standard an Strassen matrix multiplication, and LU factorization. We demonstrate that the complexity of each of these operations over the BDD representation is no greater than that over any standard representation. Further, we demonstrate that complete pivoting is no more difficult over these matrices than partial pivoting. Finally, we consider an example, the Walsh Spectrum of a Boolean function.

On the Parallel Complexity of Gaussian Elimination with Pivoting

Consider the Gaussian elimination algorithm with the well-known partial pivoting strategy for improving numerical stability (GEPP). Vavasis proved that the problem of determining the pivot sequence used by GEPP is log space-complete forP, and thus inherently sequential. AssumingP≠C, we prove here that either the latter problem cannot be solved in parallel timeO(N1/2−ε) or all the problems inPadmit polynomial speedup. HereNis the order of the input matrix andεis any positive constant. This strengthens the P-completeness result mentioned above. We conjecture that the result proved in this paper holds for the stronger boundO(N1−ε) as well, and provide supporting evidence for the conjecture. Note that this is equivalent to asserting the asymptotic optimality of the naive parallel algorithm for GEPP (moduloP≠NC).

IMPLICIT GAUSS-JORDAN SCHEME FOR THE SOLUTION OF LINEAR SYSTEMS

In this paper we discuss the results of the parallel implementation of a new scheme for the solution of large and dense linear systems: an implicit version of the well known Gauss-Jordan method. We present a complete computational analysis of the method and we explain the superiority of implicit schemes on the basis of the number of accesses to the problem data. A new pivoting strategy for the implicit methods is presented. We prove that this strategy gives better results than the classical partial pivoting used with Gaussian Elimination.