Large-scale optimization algorithms frequently require sparse Hessian matrices that are not readily available. Existing methods for approximating large sparse Hessian matrices either do not impose sparsity or are computationally prohibitive. To try and overcome these limitations, we propose a novel approach that seeks to satisfy as many componentwise secant equations as necessary to define each row of the Hessian matrix. A naive application of this approach is too expensive for Hessian matrices that have some relatively dense rows but, by carefully taking into account the symmetry and connectivity of the Hessian matrix, we are able devise an approximation algorithm that is fast and efficient with scope for parallelism. Example sparse Hessian matrices from the CUTEst test collection for optimization illustrate the effectiveness and robustness of our proposed method.

The emergence of low precision floating-point arithmetic in computer hardware has led to a resurgence of interest in the use of mixed precision numerical linear algebra. For linear systems of equations, there has been renewed enthusiasm for mixed precision variants of iterative refinement. We consider the iterative solution of large sparse systems using incomplete factorization preconditioners. The focus is on the robust computation of such preconditioners in half precision arithmetic and employing them to solve symmetric positive definite systems to higher precision accuracy; however, the proposed ideas can be applied more generally. Even for well-conditioned problems, incomplete factorizations can break down when small entries occur on the diagonal during the factorization. When using half precision arithmetic, overflows are an additional possible source of breakdown. We examine how breakdowns can be avoided and implement our strategies within new half precision Fortran sparse incomplete Cholesky factorization software. Results are reported for a range of problems from practical applications. These demonstrate that, even for highly ill-conditioned problems, half precision preconditioners can potentially replace double precision preconditioners, although unsurprisingly this may be at the cost of additional iterations of a Krylov solver.

Large-scale overdetermined linear least squares problems arise in many practical applications. One popular solution method is based on the backward stable QR factorization of the system matrix A . This article focuses on sparse-dense least squares problems in which A is sparse except from a small number of rows that are considered dense. For large-scale problems, the direct application of a QR solver either fails because of insufficient memory or is unacceptably slow. We study several solution approaches based on using a sparse QR solver without modification, focussing on the case that the sparse part of A is rank deficient. We discuss partial matrix stretching and regularization and propose extending the augmented system formulation with iterative refinement for sparse problems to sparse-dense problems, optionally incorporating multi-precision arithmetic. In summary, our computational study shows that, before applying a black-box QR factorization, a check should be made for rows that are classified as dense and, if such rows are identified, then A should be split into sparse and dense blocks; a number of ways to use a black-box QR factorization to exploit this splitting are possible, with no single method found to be the best in all cases.

We recently introduced a sparse stretching strategy for handling dense rows that can arise in large-scale linear least-squares problems and make such problems challenging to solve. Sparse stretching is designed to limit the amount of fill within the stretched normal matrix and hence within the subsequent Cholesky factorization. While preliminary results demonstrated that sparse stretching performs significantly better than standard stretching, it has a number of limitations. In this article, we discuss and illustrate these limitations and propose new strategies that are designed to overcome them. Numerical experiments on problems arising from practical applications are used to demonstrate the effectiveness of these new ideas. We consider both direct and preconditioned iterative solvers.

In recent years, a variety of preconditioners have been proposed for use in solving large sparse linear least-squares problems. These include simple diagonal preconditioning, preconditioners based on incomplete factorizations, and stationary inner iterations used with Krylov subspace methods. In this study, we briefly review preconditioners for which software has been made available, then present a numerical evaluation of them using performance profiles and a large set of problems arising from practical applications. Comparisons are made with state-of-the-art sparse direct methods.

Sparse symmetric indefinite problems arise in a large number of important application areas; they are often solved through the use of an LDL T factorization via a sparse direct solver. While for many problems prescaling the system matrix A is sufficient to maintain stability of the factorization, for a small but important fraction of problems numerical pivoting is required. Pivoting often incurs a significant overhead, and consequently, a number of techniques have been proposed to try and limit the need for pivoting. In particular, numerically aware ordering algorithms may be used, that is, orderings that depend not only on the sparsity pattern of A but also on the values of its (scaled) entries. Current approaches identify large entries of A and symmetrically permute them onto the subdiagonal, where they can be used as part of a 2 × 2 pivot. This is numerically effective, but the fill in the factor L and hence the runtime of the factorization and subsequent triangular solves may be significantly increased over a standard ordering if no pivoting is required. We present a new algorithm that combines a matching-based approach with a numerically aware nested dissection ordering. Numerical comparisons with current approaches for some tough symmetric indefinite problems are given.

In recent years, performance profiles have become a popular and widely used tool for benchmarking and evaluating the performance of several solvers when run on a large test set. Here we use data from a real application as well as a simple artificial example to illustrate that caution should be exercised when trying to interpret performance profiles to assess the relative performance of the solvers.

In recent years, there has been considerable interest in the potential for graphics processing units (GPUs) to speed up the performance of sparse direct linear solvers. Efforts have focused on symmetric positive-definite systems for which no pivoting is required, while little progress has been reported for the much harder indefinite case. We address this challenge by designing and developing a sparse symmetric indefinite solver SSIDS. This new library-quality LDL T factorization is designed for use on GPU architectures and incorporates threshold partial pivoting within a multifrontal approach. Both the factorize and the solve phases are performed using the GPU. Another important feature is that the solver produces bit-compatible results. Numerical results for indefinite problems arising from a range of practical applications demonstrate that, for large problems, SSIDS achieves performance improvements of up to a factor of 4.6 × compared with a state-of-the-art multifrontal solver on a multicore CPU.

This article focuses on the design and development of a new robust and efficient general-purpose incomplete Cholesky factorization package HSL_MI28, which is available within the HSL mathematical software library. It implements a limited memory approach that exploits ideas from the positive semidefinite Tismenetsky-Kaporin modification scheme and, through the incorporation of intermediate memory, is a generalization of the widely used ICFS algorithm of Lin and Moré. Both the density of the incomplete factor and the amount of memory used in its computation are under the user's control. The performance of HSL_MI28 is demonstrated using extensive numerical experiments involving a large set of test problems arising from a wide range of real-world applications. The numerical experiments are used to isolate the effects of scaling, ordering, and dropping strategies so as to assess their usefulness in the development of robust algebraic incomplete factorization preconditioners and to select default settings for HSL_MI28. They also illustrate the significant advantage of employing a modest amount of intermediate memory. Furthermore, the results demonstrate that, with limited memory, high-quality yet sparse general-purpose preconditioners are obtained. Comparisons are made with ICFS, with a level-based incomplete factorization code and, finally, with a state-of-the-art direct solver.

The performance of a sparse direct solver is dependent upon the pivot sequence that is chosen before the factorization begins. In the case of symmetric indefinite systems, it may be necessary to modify this sequence during the factorization to ensure numerical stability. These modifications can have serious consequences in terms of time as well as the memory and flops required for the factorization and subsequent solves. This study focuses on hard-to-solve sparse symmetric indefinite problems for which standard threshold partial pivoting leads to significant modifications. We perform a detailed review of pivoting strategies that are aimed at reducing the modifications without compromising numerical stability. Extensive numerical experiments are performed on a set of tough problems arising from practical applications. Based on our findings, we make recommendations on which strategy to use and, in particular, a matching-based approach is recommended for numerically challenging problems.

At the heart of a frontal or multifrontal solver for the solution of sparse symmetric sets of linear equations, there is the need to partially factorize dense matrices (the frontal matrices) and to be able to use their factorizations in subsequent forward and backward substitutions. For a large problem, packing (holding only the lower or upper triangular part) is important to save memory. It has long been recognized that blocking is the key to efficiency and this has become particularly relevant on modern hardware. For stability in the indefinite case, the use of interchanges and 2 × 2 pivots as well as 1 × 1 pivots is equally well established. In this article, the challenge of using these three ideas (packing, blocking, and pivoting) together is addressed to achieve stable factorizations of large real-world symmetric indefinite problems with good execution speed. The ideas are not restricted to frontal and multifrontal solvers and are applicable whenever partial or complete factorizations of dense symmetric indefinite matrices are needed.

On many current and emerging computing architectures, single-precision calculations are at least twice as fast as double-precision calculations. In addition, the use of single precision may reduce pressure on memory bandwidth. The penalty for using single precision for the solution of linear systems is a potential loss of accuracy in the computed solutions. For sparse linear systems, the use of mixed precision in which double-precision iterative methods are preconditioned by a single-precision factorization can enable the recovery of high-precision solutions more quickly and use less memory than a sparse direct solver run using double-precision arithmetic. In this article, we consider the use of single precision within direct solvers for sparse symmetric linear systems, exploiting both the reduction in memory requirements and the performance gains. We develop a practical algorithm to apply a mixed-precision approach and suggest parameters and techniques to minimize the number of solves required by the iterative recovery process. These experiments provide the basis for our new code HSL_MA79—a fast, robust, mixed-precision sparse symmetric solver that is included in the mathematical software library HSL. Numerical results for a wide range of problems from practical applications are presented.

Out-of-core sparse direct solvers reduce the amount of main memory needed to factorize and solve large sparse linear systems of equations by holding the matrix data, the computed factors, and some of the work arrays in files on disk. The efficiency of the factorization and solution phases is dependent upon the number of entries in the factors. For a given pivot sequence, the level of fill in the factors beyond that predicted on the basis of the sparsity pattern alone depends on the number of pivots that are delayed (i.e., the number of pivots that are used later than expected because of numerical stability considerations). Our aim is to limit the number of delayed pivots, while maintaining robustness and accuracy. In this article, we consider a new out-of-core multifrontal solver HSL_MA78 from the HSL mathematical software library that is designed to solve the unsymmetric sparse linear systems that arise from finite element applications. We consider how equilibration can be built into the solver without requiring the system matrix to be held in main memory. We also examine the effects of different pivoting strategies, including threshold partial pivoting, threshold rook pivoting, and static pivoting. Numerical experiments on problems arising from a range of practical applications illustrate the importance of scaling and show that, in some cases, rook pivoting can be more efficient than partial pivoting in terms of both the factorization time and the sparsity of the computed factors.

Fortran_Virtual_Memory is a Fortran 95 package that provides facilities for reading from and writing to direct-access files. A buffer is used to avoid actual input/output operations whenever possible. The data may be spread over many files and for very large data sets these may be held on more than one device. We describe the design of Fortran_Virtual_Memory and comment on its use within an out-of-core sparse direct solver.

Direct methods for solving large sparse linear systems of equations are popular because of their generality and robustness. Their main weakness is that the memory they require usually increases rapidly with problem size. We discuss the design and development of the first release of a new symmetric direct solver that aims to circumvent this limitation by allowing the system matrix, intermediate data, and the matrix factors to be stored externally. The code, which is written in Fortran and called HSL_MA77, implements a multifrontal algorithm. The first release is for positive-definite systems and performs a Cholesky factorization. Special attention is paid to the use of efficient dense linear algebra kernel codes that handle the full-matrix operations on the frontal matrix and to the input/output operations. The input/output operations are performed using a separate package that provides a virtual-memory system and allows the data to be spread over many files; for very large problems these may be held on more than one device. Numerical results are presented for a collection of 30 large real-world problems, all of which were solved successfully.

In recent years a number of solvers for the direct solution of large sparse symmetric linear systems of equations have been developed. These include solvers that are designed for the solution of positive definite systems as well as those that are principally intended for solving indefinite problems. In this study, we use performance profiles as a tool for evaluating and comparing the performance of serial sparse direct solvers on an extensive set of symmetric test problems taken from a range of practical applications.

We recently carried out an extensive comparison of the performance of state-of-the-art sparse direct solvers for the numerical solution of symmetric linear systems of equations. Some of these solvers were written primarily as research codes while others have been developed for commercial use. Our experiences of using the different packages to solve a wide range of problems arising from real applications were mixed. In this paper, we highlight some of these experiences with the aim of providing advice to both software developers and users of sparse direct solvers. We discuss key features that a direct solver should offer and conclude that while performance is an essential factor to consider when choosing a code, there are other features that a user should also consider looking for that vary significantly between packages.

The need to solve large sparse linear systems of equations efficiently lies at the heart of many applications in computational science and engineering. For very large systems when using direct factorization methods of solution, it can be beneficial and sometimes necessary to use multiple processors, because of increased memory availability as well as reduced factorization time. We report on the development of a new parallel code that is designed to solve linear systems with a highly unsymmetric sparsity structure using a modest number of processors (typically up to about 16). The problem is first subdivided into a number of loosely connected subproblems and a variant of sparse Gaussian elimination is then applied to each of the subproblems in parallel. An interface problem in the variables on the boundaries of the subproblems must also be factorized. We discuss how our software is designed to achieve the goals of portability, ease of use, efficiency, and flexibility, and illustrate its performance on an SGI Origin 2000, a Cray T3E, and a 2-processor Compaq DS20, using problems arising from real applications.

In recent years, a number of new direct solvers for the solution of large sparse, symmetric linear systems of equations have been added to the mathematical software library HSL. These include solvers that are designed for the solution of positive-definite systems as well as solvers that are principally intended for solving indefinite problems. The available choice can make it difficult for users to know which solver is the most appropriate for their use. In this study, we use performance profiles as a tool for evaluating and comparing the performance of the HSL solvers on an extensive set of test problems taken from a range of practical applications.

Many applications in science and engineering give rise to large sparse linear systems of equations that need to be solved as efficiently as possible. As the size of the problems of interest increases, it can become necessary to consider exploiting multiprocessors to solve these systems. We report on the design and development of parallel frontal solvers for the numerical solution of large sparse linear systems. Three codes have been developed for the mathematical software library HSL (www.cse.clrc.ac.uk/Activity/HSL). The first is for unsymmetric finite-element problems; the second is for symmetric positive definite finite-element problems; and the third is for highly unsymmetric linear systems such as those that arise in chemical process engineering. In each case, the problem is subdivided into a small number of loosely connected subproblems and a frontal method is then applied to each of the subproblems in parallel. We discuss how our software is designed to achieve the goals of portability, ease of use, efficiency, and flexibility, and illustrate the performance using problems arising from real applications.

Hager recently introduced down and up exchange methods for reducing the profile of a sparse matrix with a symmetric sparsity pattern. The methods are particularly useful for refining orderings that have been obtained using a standard profile reduction algorithm, such as the Sloan method. The running times for the exchange algorithms reported by Hager suggested their cost could be prohibitive for practical applications. We examine how to implement the exchange algorithms efficiently. For a range of real test problems, it is shown that the cost of running our new implementation does not add a prohibitive overhead to the cost of the original reordering.

We describe the design, implementation, and performance of a frontal code for the solution of large sparse symmetric systems of linear finite-element equations. The code is intended primarily for positive-definite systems, since numerical pivoting is not performed. The resulting software package, MA62, will be included in the Harwell Subroutine Library. We illustrate the performance of our new code on a range of problems arising from real engineering and industrial applications. The performance of the code is compared with that of the Harwell Subroutine Library general frontal solver MA42 and with other positive-definite codes from the Harwell Subroutine Library.

We describe the design, implementation, and performance of a frontal code for the solution of large, sparse, unsymmetric systems of linear equations. The resulting software package, MA42, is included in Release 11 of the Harwell Subroutine Library and is intended to supersede the earlier MA32 package. We discuss in detail the extensive use of higher-level BLAS kernels within MA42 and illustrate the performance on a range of practical problems on a CRAY Y-MP, an IBM 3090, and an IBM RISC System/6000. We examine extending the frontal solution scheme to use multiple fronts to allow MA42 to be run in parallel. We indicate some directions for future development.

Arnoldi methods can be more effective than subspace iteration methods for computing the dominant eigenvalues of a large, sparse, real, unsymmetric matrix. A code, EB12 , for the sparse, unsymmetric eigenvalue problem based on a subspace iteration algorithm, optionally combined with Chebychev acceleration, has recently been described by Duff and Scott and is included in the Harwell Subroutine Library. In this article we consider variants of the method of Arnoldi and discuss the design and development of a code to implement these methods. The new code, which is called EB13 , offers the user the choice of a basic Arnoldi algorithm, an Arnoldi algorithm with Chebychev acceleration, and a Chebychev preconditioned Arnoldi algorithm. Each method is available in blocked and unblocked form. The code may be used to compute either the rightmost eigenvalues, the eigenvalues of largest absolute value, or the eigenvalues of largest imaginary part. The performance of each option in the EB13 package is compared with that of subspace iteration on a range of test problems, and on the basis of the results, advice is offered to the user on the appropriate choice of method. — Author's Abstract

This paper discusses the design and development of a code to calculate the eigenvalues of a large sparse real unsymmetric matrix that are the rightmost, leftmost, or are of the largest modulus. A subspace iteration algorithm is used to compute a sequence of sets of vectors that converge to an orthonormal basis for the invariant subspace corresponding to the required eigenvalues. This algorithm is combined with Chebychev acceleration if the rightmost or leftmost eigenvalues are sought, or if the eigenvalues of largest modulus are known to be the rightmost or leftmost eigenvalues. An option exists for computing the corresponding eigenvectors. The code does not need the matrix explicitly since it only requires the user to multiply sets of vectors by the matrix. Sophisticated and novel iteration controls, stopping criteria, and restart facilities are provided. The code is shown to be efficient and competitive on a range of test problems.