Pivoting Research Articles

A partial pivoting strategy for sparse symmetric matrix decomposition

It is well known that the partial pivoting strategy by Bunch and Kaufman is very effective for factoring dense symmetric indefinite matrices using the diagonal pivoting method. In this paper, we study a threshold version of the strategy for sparse symmetric matrix decomposition. The use of this scheme is explored in the multifrontal method of Duff and Reid for sparse indefinite systems. Experimental results show that it is at least as effective as the existing pivoting strategy used in the current multifrontal implementation.

The effect of partial pivoting in sparse ordinary differential equation solvers

This paper considers the effect of partial pivoting in automatic ordinary differential equation solvers that utilize sparse matrix techniques. Two solvers are considered, the well-known LSODES solver and a derivate LSOD28. LSODES uses the Yale Sparse Matrix Package which does not perform partial pivoting. LSOD28 uses the MA28 Sparse Matrix Package which does perform partial pivoting. Results are presented for a benchmark problem that contains several features typically present in realistic problems. The results demonstrate that both solvers perform satisfactorily. At the same time, they illustrate that the lack of partial pivoting does not necessarily degrade the efficiency or the reliability of a solver such as LSODES for such problems. They support the argument that partial pivoting is not necessarily required in adaptive sparse solvers to solve complex problems accurately and efficiently.

Condition Number Estimators in a Sparse Matrix Software

The stability of the computational process in the solution of systems of linear algebraic equations $Ax = b$ depends on the condition number of matrix A. Reliable and efficient algorithms for calculating estimates of the condition number of a matrix are given in An estimate for the condition number of a matrix, SIAM J. Numer. Anal., 16(1979), pp. 368–375. The application of these algorithms in a sparse matrix software, package Y12M [Y12M-solution of large and sparse systems of linear algebraic equations, Lecture Notes in Computer Science, Vol. 121, Springer, Berlin, 1981], is discussed. Three algorithms have been implemented in package Y12M and tested on a very large set of problems. The pivotal strategies for sparse matrices, which are used instead of the partial pivoting for dense matrices, usually provide reliable estimates for the condition number when the value of the stability factor u is small, say $u \in [4,16]$ ($u = 1$ for the partial pivoting). It is shown that for the subroutines of package Y12M this is true even if the stability factor u is rather large. This phenomenon is explained by a careful analysis of the main pivotal strategy in package Y12M. Very often a special nonnegative parameter, a drop-tolerance, is used during the factorization of matrix A so that each element, which in the course of the calculations becomes smaller in absolute value than the drop-tolerance, is considered as a zero element. Both storage and computing time may be saved in this way for some classes of matrices. The influence of the use of a large drop-tolerance on the reliability of the condition number estimators is discussed. Some conclusions, concerning the three condition number estimators, are given.

Column LU Factorization with Pivoting on a Message-Passing Multiprocessor

A column-oriented algorithm is presented for LU factorization with partial pivoting. Two different mappings of columns to processors are considered. Forward and backsubstitution algorithms to use the factorization for solving linear systems are developed. Timing data, including processor utilization and load balancing, are provided by a hypercube simulator.

Pipe Network Analysis by Partial Pivoting Method

A simple and efficient, computer‐oriented method of analyzing hydraulic networks is presented, in which the nonlinear energy equations are linearized to yield, together with the continuity equations, a set of equations, [A][Q]=[C], that is solved by Guassian elimination process refined with partial pivoting. The method is reliable, and is shown to compare favorably with others.

Inverse laminar boundary‐layer problems with assigned shear. The mechul function revisited

AbstractAn inverse method is presented which accurately determines the pressure distribution for assigned wall shear in a two‐dimensional, laminar, incompressible boundary layer. The method reformulates the mechul function scheme of Cebeci and Keller to produce a stable solution in the marching direction and to increase accuracy in the normal direction. In the reformulation a modified pressure gradient parameter variation in the normal direction is used in conjunction with three‐point backward differences for streamwise derivatives and fourth‐order accurate splines for normal derivatives. The resulting spline‐finite difference equations are solved by Newton‐Raphson iteration together with partial pivoting. Numerical solutions are presented for self‐similar and non self‐similar flows and compared with published results.

Three mysteries of Gaussian elimination

If numerical analysts understand anything, surely it must be Gaussian elimination. This is the oldest and truest of numerical algorithms. To be precise, I am speaking of Gaussian elimination with partial pivoting, the universal method for solving a dense, unstructured n X n linear system of equations Ax = b on a serial computer. This algorithm has been so successful that to many of us, Gaussian elimination and Ax = b are more or less synonymous. The chapter headings in the book by Golub and Van Loan [3] are typical -- along with "Orthogonalization and Least Squares Methods," "The Symetric Eigenvalue Problem," and the rest, one finds "Gaussian Elimination," not "Linear Systems of Equations."

Numerical Solution of Two-Dimensional Axisymmetric Free Convection within Isotope Separating Thermal Diffusion Column

A numerical method for integrating the equations of change without the Boussinesque approximation is presented for axisymmetric free convection within thermal diffusion column. A finite difference analogue of equations of change in terms of primitive variables (density, radial and axial velocities, and temperature) is solved by a modified Newton method. Errors in a finite difference approximation of derivatives for temperature and density are reduced by introduction of a perturbation or deviation from an analytically obtainable base temperature distribution. Roundoff errors in a Gaussian elimination procedure are also reduced with the adoption of “scaling” and a partial pivoting strategy of the Jacobian matrix. A relation for total mass conservation is used to prevent the density from converging to a trivial solution, i.e. zero. The method has been applied to a total reflux flow of Ar gas within the 1,000 mm height thermal diffusion column with an inner hot radius of 0.2 mm and an outer cold radius of 5 mm, between which the temperature differs by 50∼400K. The “end effect” appears within 15 mm portion from the top and bottom plates.

An Implementation of Gaussian Elimination with Partial Pivoting for Sparse Systems

In this paper, we consider the problem of solving a sparse nonsingular system of linear equations. We show that the structures of the triangular matrices obtained in the $LU$-decomposition of a sparse nonsingular matrix A using Gaussian elimination with partial pivoting are contained in those of the Cholesky factors of $A^T A$, provided that the diagonal elements of A are nonzero. Based on this result, a method for solving sparse linear systems is then described. The main advantage of this method is that the numerical computation can be carried out using a static data structure. Numerical experiments comparing this method with other implementations of Gaussian elimination for solving sparse linear systems are presented and the results indicate that the method proposed in this paper is quite competitive with other approaches.

A note on estimating the error in Gaussian elimination without pivoting

This article deals with the problem of estimating the error in the computed solution to a system of equations when that solution is obtained by using Gaussian elimination without pivoting. The corresponding problem, where either partial or complete pivoting is used, has received considerable attention, and efficient and reliable methods have been developed. However, in the context of solving large sparse systems, it is often very attractive to apply Gaussian elimination without pivoting, even though it cannot be guaranteed a-priori that the computation is numerically stable. When this is done, it is important to be able to determine when serious numerical errors have occurred, and to be able to estimate the error in the computed solution. In this paper a method for achieving this goal is described. Results of a large number of numerical experiments suggest that the method is both inexpensive and reliable.

Numerical solution of two-dimensional axisymmetric free convection within isotope separating thermal diffusion column.

A numerical method for integrating the equations of change without the Boussinesque approximation is presented for axisymmetric free convection within thermal diffusion column. A finite difference analogue of equations of change in terms of primitive variables (density, radial and axial velocities, and temperature) is solved by a modified Newton method. Errors in a finite difference approximation of derivatives for temperature and density are reduced by introduction of a perturbation or deviation from an analytically obtainable base temperature distribution. Roundoff errors in a Gaussian elimination procedure are also reduced with the adoption of “scaling” and a partial pivoting strategy of the Jacobian matrix. A relation for total mass conservation is used to prevent the density from converging to a trivial solution, i.e. zero. The method has been applied to a total reflux flow of Ar gas within the 1,000 mm height thermal diffusion column with an inner hot radius of 0.2 mm and an outer cold radius of 5 m...

Practical Use of Some Krylov Subspace Methods for Solving Indefinite and Nonsymmetric Linear Systems

The main purpose of this paper is to develop stable versions of some Krylov subspace methods for solving linear systems of equations $Ax = b$. As in the case of Paige and Saunders’s SYMMLQ [SIAM J. Numer. Anal., 12 (1975), pp. 617–624], our algorithms are based on stable factorizations of the banded Hessenberg matrix representing the restriction of the linear application A to a Krylov subspace. We will show how an algorithm similar to the SYMMLQ can be derived for nonsymmetric problems and we will describe a more economical algorithm based upon the $LU$ factorization with partial pivoting. In the particular case where A is symmetric indefinite the new algorithm is theoretically equivalent to SYMMLQ but slightly more economical. As a consequence, an advantage of the new approach is that nonsymmetric or symmetric indefinite or both nonsymmetric and indefinite systems of linear equations can be handled by a single algorithm.

Collocation for Singular Perturbation Problems I: First Order Systems with Constant Coefficients

The application of collocation methods for the numerical solution of singularly perturbed ordinary differential equations is investigated. Collocation at Gauss, Radau and Lobatto points is considered, for both initial and boundary value problems for first order systems with constant coefficients. Particular attention is paid to symmetric schemes for boundary value problems; these problems may have boundary layers at both interval ends. .br Our analysis shows that certain collocation schemes, in particular those based on Gauss or Lobatto points, do perform very well on such problems, provided that a fine mesh with steps proportional to the layers'' width is used in the layers only, and a coarse mesh, just fine enough to resolve the solution of the reduced problem, is used in between. Ways to construct appropriate layer meshes are proposed. Of all methods considered, the Lobatto schemes appear to be the most promising class of methods, as they essentially retain their usual superconvergence power for the smooth, reduced solution, whereas Gauss-Legendre schemes do not. .br We also investigate the conditioning of the linear systems of equations arizing in the discretization of the boundary value problem. For a row equilibrated version of the discretized system we obtain a pleasantly small bound on the maximum norm condition number, which indicates that these systems can be solved safely by Gaussian elimination with scaled partial pivoting.

Effect of Equilibration on Residual Size for Partial Pivoting

It is shown that column pivoting with row equilibration satisfies the same type of error bound as does row pivoting without scaling, and that row pivoting with column equilibration satisfies the same type of bound as does column pivoting without scaling. An interesting consequence for row pivoting is that column equilibration is sufficient to ensure that the norm of the residual is reasonably small.

Computations of coefficients in the polynomials of Pade´approximations by solving systems of linear equations

We explore reliability, stability and accuracy of determining the polynomials which define the Pade´approximation to a given function h(x) by solving a system of linear equations to get the coefficients in the denominator polynomial B n(x). The coefficients in the numerator polynomial A m(x) follow directly from those for B n(x). Our approach is in the main heuristic. For the numerics we use the models e −x1n(1 +x), (1 +x) ± 1/2 and the exponential integral, each with m=n. The system of equations, with matrix of Toeplitz type, was solved by Gaussian elimination (Crout algorithm) with equilibration and partial pivoting. For each model, the maximum number of incorrect figures in the coefficients is of the order n at least, thus indicating that the matrix becomes ill conditioned as n increases. Let δ n(x)andω n(x) be the errors in A n(x) and B n(x) respectively, due to errors in the coefficients of B n(x). For x fixed, δ n(x) and ω n(x) and the corresponding relative errors increase as n increases. However, for a considerable range on n, the relative errors in A n(x)B n(x) are virtually nil. This has the following theoretical explanation. Now B n(x)h(x) −A m(m) = 0 (x m+n+ 1). It can be shown that ω n(x)h(x) − δ m(x) = 0(x m+ 1). In this sense both A m(x)B n(x)andδ m(x)ω n(x) are approximations to h(x). Thus if the difference of these two approximations and ω n(x)B n(x), the relative error in B n(x), are sufficiently small, then the relative error in A m(x)/B n(x) is of no consequence.

On factoring a class of complex symmetric matrices without pivoting

Let A = B + i C \mathcal {A} = \mathcal {B} + i\mathcal {C} be a complex, symmetric n × n n \times n matrix with B \mathcal {B} and C \mathcal {C} each real, symmetric and positive definite. We show that the LINPACK diagonal pivoting decomposition U − 1 A ( U − 1 ) T = D {\mathcal {U}^{ - 1}}\mathcal {A}{({\mathcal {U}^{ - 1}})^T} = \mathcal {D} proceeds without the necessity for pivoting. In particular, when B \mathcal {B} and C \mathcal {C} are band matrices, bandwidth is preserved.

Elastic-Plastic Tension-Torsion in a Circular Bar of Rate-Sensitive Material

A theoretical and numerical investigation is made of the behavior of a circular bar of elastic-viscoplastic material subjected to either proportional or nonproportional straining in tension and torsion. The tensile properties of the material are assumed to be adequately represented by the empirical formula of Bodner and Symonds, which is based on experimental results for mild steel. The analysis employs Perzyna’s generalization of von Mises’ yield criterion, and the associated plastic flow rule. Elastic compressibility is taken into account. The stresses, strains, and loading paths are determined for four prescribed straining paths; the results were obtained by solving numerically, using Crout’s method with partial pivoting, a system of six simultaneous quasi-linear partial differential equations of hyperbolic type. Closed-form analytical solutions are also obtained for the stresses, load, and torque corresponding to fully plastic conditions; interaction curves are plotted relating load and torque under these conditions, for quasi-static straining and for a range of finite effective strain rates.

Algorithms for sparse Gaussian elimination with partial pivoting

Algorithms for sparse Gaussian elimination with partial pivoting

Two algorithms for treating ax = λbx

The first algorithm is a generalization of the iterative method of successive coordinate overrelaxation. Instead of computing a single eigenvalue and its corresponding eigenvector at a time, p vectors are iterated simultaneously in such a way that they converge ultimately towards the eigenvectors corresponding to the p smallest eigenvalues. This method of simultaneous group coordinate overrelaxation has the advantage of taking into account the sparseness of the matrices A and B to full extent, to operate on the originally given matrices and to require a relatively small working storage space. The second algorithm is a variant of the well-known bisection method. It is shown that the inertia theorem for quadratic forms can be applied to localize the eigenvalues. Hence, an appropriate congruence transformation of the matrix A − β B , which is assumed now to be banded, allows one to determine the number of eigenvalues that are smaller than μ. In order to increase the numerical stability of the reduction process, a partial pivoting is applied consisting of auxiliary congruence transformations in such a way that the band width is only doubled. This approach to the bisection method is conceptually much simpler - the actual algorithm preserves the symmetry and requires a smaller working storage space for the computation of the eigenvectors by inverse vector iteration in comparison with the well-known realizations of the bisection method.

A note on partial pivoting and Gaussian elimination

The results of Wilkinson on the stability of Gaussian elimination with column pivoting are generalized.