EFFICIENT SOLUTION OF LINEAR SPARSE SYSTEMS

Abstract

The present thesis is devoted for finding the efficient solution of sparse systems of linear equations. One of the most frequently encountered problems in scientific computing is the manipulation and solution of large sparse unstructured (or structured) systems of linear equations Ax = b , where A is a non-singular square matrix of order n and b is a given column vector of order 17. The term sparse means most of the entries in the matrix A are zeros. The application of finite difference, finite element and domain decomposition methods etc. to approximate the solution of partial differential equations produces a sparse system of linear equations. Much research in the past thirty years has been targeted for solving these systems and several sophisticated techniques and codes are made available for solving them. Computational techniques for solving a sparse system of linear equations may be readily classified into either direct or iterative methods. With few exceptions all the methods exploit the sparsity of the coefficient matrix A . In the former category a solution is guaranteed on the completion of a fixed amount of arithmetic work, whilst the iterative methods generally involve a repetitive sequence of simple matrix-vector operations in which a guess vector is successfully improved until the solution is obtained to a specified accuracy. In the iterative solution a sequence of approximations is generated iteratively which (hopefully) converges to a solution of the linear system. However, for certain classes of problems, use of an appropriate iterative method can yield an approximate solution significantly faster than direct methods. Unlike direct methods, iterative methods tend to be more special-purpose and it is well known that there are no effective iterative algorithms for unstructured (arbitrary) sparse system of linear equations. However, the iterative methods may not converge or their convergence rate can be very slow which could suggest that direct methods should be preferred. Several direct methods are available in literature for solving a sparse system of linear equations. The sparse LU-decomposition method is one of the most popular directed methods for solving Ax = b with different b's. An important concept in sparse LU-decomposition is the determination of the sparsity pattern of the factors L and U . Moreover, the determination of fill- ins of a sparse matrix is a central problem in the solution of Ax =b . The appearance of the fill-ins ir. LU-decomposition is not the only trouble. The major problem is that the positions in which appear are not known in advance. This problem is discussed as one of the objectives of the thesis. The Krylov subspaces are some of the most useful and powerful tools for computing approximate solutions in different Linear Algebra problems. Some of the most efficient iterative methods are those of Krylov subspace type. Iterative methods in this class of algorithms construct an approximate solution to the problem Ax b in the Krylov subspace. The solution to a linear system of equations lies in a Krylov space whose dimension is the degree of the minimal polynomial of A . Therefore, if the minimal polynomial of A has low degree then the space in which a Krylov method searches for the solution can be small. In this case a Krylov method has opportunity to converge fast. Although the additional storage for iterative methods is small, in all but the simplest cases they do not work: convergence is slower or non-existent for problems coming from non-trivial applications. To have any chance of beating the direct methods, the iterative methods must converge rapidly, and this naturally leads to the search for good preconditioners for A . A simple way to do this is to replace the original problem Ax = b with a new, easier problem with same solution. A preconditioner is a non-singular matrix M whose inverse is assumed to be a close approximation to the inverse of A and it should be chosen such that some convenient numerical procedure is available for solving the preconditioning step Mz = v. Several preconditioning techniques have been proposed in literature for solving sparse linear systems of equations. The incomplete factorization-based preconditioners are some of the most efficient preconditioning techniques for solving Ax = b . The TW(l ) and ILUT( k0, z) are such preconditioners. The parameter 1 denotes the level of fill-in, which is to be allowed in the factorization. k0 is an integer such that each row of L and U will have a maximum of k0 elements in addition to their original number of non-zero elements, z- is some tolerance used for dropping entries in L and U. The incomplete LU( l) preconditioner uses an incomplete LU-decomposition of the matrix A as a preconditioner, say M = LU. To be effective, the ILU( ) oreconditioner requires the determination of an efficient technique for finding the sparsity pattern of L ii and U . In this thesis, a good method for determining the sparsity pattern of the factors L and U before incomplete LU-decomposition usually starts, is presented. The whole range of the subject of work presented in the thesis has been divided into eight Chapters. Chapterwise summary of the thesis is as follows. Chapter 1 The Chapter 1 is of introductory nature. In this Chapter, two different types of methods-direct methods and iterative methods to solve a sparse system of linear equations are presented. The details of relevant numerical methods and a brief account of the related studies done by previous researchers are given. The presentation of storage scheme for saving input data of a sparse matrix is also discussed. A brief introduction of the work done in other Chapters is also given. The libraries of IIT, Roorkee; IIT, New Delhi; Indian Statistical Institute, New Delhi; University of Delhi and personal efforts by requesting authors of the papers have been the provision sources for the study material. Chapter 2 The Chapter 2 is devoted for finding the direct solution of a sparse system of linear equations by using Cramer's rule based on digraph approach. Cramer's rule is one of the well-known direct methods for solving Ax = b by evaluating of a numerous number of determinants. Cramer's rule produces the solution directly without creating fill-ins encountered in other direct methods. While the Cramer's rule is an important tool in many theoretical studies of the solution of linear systems of equations, it is practically never used in practice; especially when the systems solved are large. This is due to high computational cost of this algorithm. Indeed, its cost is 0(n!) , while the normally used Gaussian's elimination can be performed in 0(n3 simple arithmetic operations. Of course, these costs are valid for dense matrices. In the case of sparse matrices many elements are zeros and this fact may lead to a dramatical reduction of the computational cost. Therefore, the idea of using the Cramer's rule for sparse matrices is interesting. We decided to concentrate on an area in linear systems that was paid a little attention in recent years and that depends heavily on the theory of digraphs and theory of trees. iii Digraph approach is one of the methods based on Cramer's rule for so lving Ax=b. Applications of digraph approach in the analysis of a linear system such as networks give expressions for the determinant of A in terms of labeled subdigraphs in digraphs. Thus, Cramer's rule can be reformulated in terms of these labeled subgraphs in digraphs. In the present Chapter, we propose an algorithm that depends on the digraph approach to solve a sparse system of linear equations. It is based on the Cramer's rule. The main feature of this algorithm is that the solution of a sparse system of linear equations can be expressed exactly if all the non-zero entries including right-hand side are integers and if all products do not exceed the size of the largest integer that can be represented in the arithmetic of the computer used. The implementation of the proposed algorithm is experimented by applying on five test problems. It is shown that efficiency with which a sparse system of linear equations can be analyzed by a digital computer using proposed digraph approach as a tool depends largely upon the efficiency with which 1-factors and 1-factorial connections are generated. Finally, the usefulness of digraph approach is discussed. The contents of this Chapter have been published in the form of a research paper in the " International Journal of Computer and Mathematics, Vol. 74, pp. 151-158, 2000 ". Chapter 3 The Chapter 3 is concentrated on the direct solution of a sparse system of linear equations by Cramer's rule based on numerical structure approach. Numerical structure technique is one of the methods based on Cramer's rule for solving Ax=b. Applications of the numerical structure in the theory of machines (such as constructing gearboxes, etc.) give expressions for the determinant of A in terms of labeled paths in trees. Thus, Cramer's rule can be reformulated in terms of these labeled paths in the trees. In Chapter 3, two algorithms in terms of tree theory are described for solving a sparse linear system of equations. The proposed algorithms depend on the numerical structure approach in which only the non-zero elements are used. The solution of the sparse system of linear equations is obtained by applying the Cramer's rule. Cramer's rule produces the solution directly without creating fill-ins encountered in other direct methods. Moreover, the solution can be expressed exactly if all the entries, including the right-hand side, are integers and if all products do not exceed the size of the largest iv integer that can be represented in the arithmetic of the computer used. We give a theorem presenting a sparsity threshold such that if the number of the non-zero entries is below this threshold the proposed algorithms will likely be more efficient than the traditionally used methods. The efficiency of the proposed algorithms for solving sparse systems of linear equations is shown by applying on seven test problems. The obtained results are also compared with that of the digraph approach described in Chapter 2. From numerical experiments, it is shown that the performance of the proposed algorithms J based on the numerical structure approach is much better than that of digraph approach. Finally, we produce a list of graph theoretical applications where the ideas given in this Chapter can be used. The application of the proposed algorithms for computing the exact permanent value of a Boolean matrix is also given. The contents of this Chapter have been published in the form of a research paper in the " Journal of Computational and Applied Mathematics, Vol. 136, pp. 1-15, 2001 ". Chapter 4 If we wish to solve Ax=b for a number of different b's, then the LU-decomposition method is more efficient than Gaussian's elimination. This is the case when performing iterative refinement and may also happen because of different analysis of a mathematical model on various inputs. Moreover, there are many applications where the entries of the inverse are useful. The sparse LU-decomposition method is the standard approach to computing entries of the inverse of a sparse matrix A. Often it is the case that a sequence of the systems of equations must be solved where the sparsity pattern is fixed but the numerical values of entry change. One example is the design problem where values of parameters must be chosen to minimize some measure of performance. We describe the best technique for determining the fill-ins, which suits these cases. The determination of fill-ins (entries that were originally zero in the matrix A can be non-zero in the factors L and U) of a sparse matrix is a central problem in the solution of sparse system of linear equations. The appearance of the fill-ins in LU-decomposition is not the only trouble. The major problem is that the positions in which fill-ins appear are not known in advance. The prediction of the fill-ins positions is the objective of the present Chapter. The basic idea of this Chapter is in determination of the pattern of non-zero entries of the LU factors of a given sparse matrix A. An efficient, effective and inexpensive new algorithm is given to determine the sparsity pattern of the factors L and U . This new algorithm is based on powers of the Boolean matrix B obtained from the coefficient matrix A such that the fill-ins in the structure of A can be known in advance. Thus, the position of fill-ins in A are determined in the best choice manner, that is, it is very effective and memory wise inexpensive. Once the sparsity pattern of L and U is obtained, the construction of the LU-decomposition is straightforward. The efficiency of the proposed sparse LU-decomposition method is shown by applying on eight test problems. Obtained results are also compared with that of the numerical structure approach of Chapter 3. The contents of this Chapter have been published in the form of a research paper in the journal " Computers and Mathematics With Applications, Vol. 43, pp. 131-155, 2002 ". Chapter 5 Several classes of Krylov subspace methods have been developed for non-symmetric problems. One of these classes is the class of restarted methods typified by the popular GMRES and GCR algorithms. Restarted method repeatedly applies k steps of GMRES method respectively to the problem Ax = b, where restarting the full minimal residual method is to reduce the build-up of an excessive number of basis vectors, which are typically very costly to the algorithm. Truncating the set of Arnoldi vectors by keeping only k vectors may lead to a loss of superlinear convergence or even to stagnation of convergence. However, since less computational work and storage are required for each iteration, the GMRES( k) method may be more efficient in practice. Although the additional storage for iterative methods typically is small (a few n -vectors), in all but the simplest cases they do not work: convergence is slower or non-existent for problems coming from non-trivial applications. To have any chance of beating the direct methods, the iterative methods must converge rapidly, and this naturally leads to the search for good preconditioners for A . Thus, in order to be effective, the GMRES( k) method must be combined with a good preconditioner. The incomplete LU-decomposition is one of the most efficient preconditioning techniques for solving Ax = b . The incomplete LU(1) preconditioner uses an incomplete LU-decomposition of the matrix A vi as a preconditioner, say M = LU . The parameter 1 denotes the level of fill-in, which is to be allowed in the factorization. To be effective, the IL U(1) preconditioner requires the determination of an efficient technique for finding the sparsity pattern of L and U . In this Chapter, a good method for determining the sparsity pattern of the factors L and U before incomplete LU-decomposition actually starts, is presented. However, the proposed method avoids the drawbacks faced earlier techniques available in literature. In particular, it does not depend upon the values of the entries of A . It is interesting, that in this method, one may choose the fill-ins to any level of factorization till complete factorization is achieved. It is assumed that, by increasing the level of fill-ins, M-1 more closely approximates A-1, and thus fewer iterations are required. In the Chapter 5, an efficient method for the construction of ILU preconditioner using powers of a Boolean matrix strategy (PBS), denoted by the symbol ILUPBS, is given. Some numerical experiments are performed which show that it is possible to considerably reduce the number of GMRES iterations if the proposed ILUPBS(1) preconditioner is used. For all tests carried out, the best value for k is found to be 10. The contents of this Chapter have been accepted for publication in the form of a research paper in the journal " Computers and Mathematics With Applications, to appear, 2002 ". Chapter 6 The GCR( k) algorithm is one of the class of the truncated methods for solving sparse non-symmetric linear systems of equations. Truncating the GCR method is to reduce the build-up of an excessive number of basis vectors, which are typically very costly to the algorithm. One will, of course, expect that the truncated GCR( k) algorithm will have a slower convergence rate to the regular procedure. However, since less computational work and storage are required for each iteration, the truncated algorithm may be more efficient in practice. The rate of convergence of the GCR(k ) algorithm can be accelerated by orthogonalization and by preconditioning. In this algorithm, preconditioning takes two forms: prior to iteration and during iteration. During the iterative process, fewer orthogonalizations means lower work requirement. The incomplete factorization-based preconditioners considered in this Chapter are the ILUPBS(l) and the vii ILUT( ko, r) preconditioners. In this Chapter, we have investigated the performance of various ILUPBS(1) and ILUT( r) preconditioned GCR( k) algorithm for different values of k, (2 k (k is number of orthogonalizations). Moreover, we make a comparison between the GCR( k) algorithm (of course for the best value of k) and the GMRES(10) solver combined with the same preconditioners. It is found that the convergence rate of ILUPBS(1)(1 1) preconditioned GCR(2) algorithm is faster than that of ILUT( ko,z- ) preconditioned GCR(2) algorithm. Our comparative simulation studied for the problems carried out here gives also the clear indication that the ILUPBS( / )( / 1) preconditioned GCR(2) solver is a good competitor of ILUPBS( / )( / 1) preconditioned GMRES(10) solver. The performance of ILUT( k0, r) preconditioned GMRES(10) method is better than that of ILUT( k0, r) preconditioned GCR(2). Finally, the performance of ILUPBS( / )( / 1) preconditioner is better than that of ILUT( k0, r) preconditioner combined with GCR(2) or GMRES(10). Chapter 7 The technique of iterative refinement for improving the computed solution to Ax=b was used on desk calculators and computers in the 1940s and has remained popular. To solve sparse linear system of equations one can try several different algorithms. One is to take a guess at the solution and iterate refining that guess until the error you make is suitably small. The method we propose is an iterative refinement based on the sparse LU-decomposition. If LU-decomposition with iterative refinement (LUIR) is used in the solution of systems of linear equations, whose matrices are dense, then the accuracy of the results will usually be greater than the accuracy obtained by the use of LU-decomposition without iterative refinement (LUDS). However, both more storage (about 100%, because a copy of matrix A is needed) and more computing time (some extra time is needed to perform the iterative process) must be used with LUIR. Normally, when the matrix A is sparse the accuracy of the solution computed by some sparse matrix techniques and LUIR will still be greater. In this Chapter, it is verified (by numerical experiments) that the use of sparse matrix techniques with LUIR may also result in a reduction of both the computing time and the viii storage requirements (this will never happen when LUIR is applied for dense matrices). The powers of the Boolean matrix strategy (PBS) is used in the efforts to achieve such a reduction. By the use of < m, (132' an attempt is made to control the sparsity of the matrix A to the level 1. The use of 1= 0,I,2,..., rrt —1 leads to an inaccurate factorization, but the accuracy lost is normally regained in the iteration process. Many test problems are given in order to compare the use of LUIR with values for 1, (m > 1 0) and the use of LUDS with m where B2in = B2m+i in the solution of systems whose matrices are large and sparse. The main conclusion is that the new iterative refinement may efficiently be used as a good option instead of those used in the packages for the solution of sparse linear systems of equations. A part of this Chapter has been accepted for presentation as a research paper in the " t h Joint 9t National Conference Of The Vijnana Of India On Applied And Industrial Mathematics" and " Ste Annual Conference Of Indi21 Society Of Information Theory And Applications", Department of Mathematics, Netaji Subhas Institute of Technology, New Delhi, India, February 22-24, 2002. Chapter 8 The Krylov subspace iterative based methods such as CGS, Bi-CGSTAB, TFQNIR and GMRES are some of the most efficient iterative methods for solving sparse non-symmetric systems of linear equations Ax = b . In order to improve the convergence of these methods, some form of preconditioning must be applied. The incomplete factorization-based preconditioners considered in this Chapter are ILUPBS( / ) and 1LUT( k0, r ). In the Chapter 8, we have presented a comparative study of the preconditioned iterative methods for solving large sparse non-symmetric linear systems of equations. The numerical implementations of these preconditioned methods are tested and compared on seven test problems. Comparisons of the effectiveness of 1LUPBS(/ ) and 1LUT( ko, z) preconditioners applied to the considered iterative methods are also shown. Numerical experiments presented herewith show that the performance of all methods combined with the ILUPBS( )(1> 0) preconditioner is better than that of them combined with the ILUT(ko, ) preconditioner. It is also found that the convergence rate of all ix methods combined with the ILUPBS(2) preconditioner is much faster than that of them combined with the ILUT( ko , r) preconditioner. Finally, numerical experiments show that there is no clear best Krylov method out of considered methods. Each of the methods is the winner in a problem and the loser in another one for the class of the test problems. A part of the present Chapter has been accepted for publication in the form of a research paper in the "Proceedings of the National Seminars on Advances in Mathematical, Statistical and Computational Methods in Sciences and Technology", Indian School of Mines, Department of Applied Mathematics, Dhanbad, India, November 29-30, 2001, To appear.