Templates for Sparse Linear Solvers

Next: Templates for Solution Up: Templates for Linear Algebra Previous: The BLAS as

# Templates for Sparse Linear Solvers

The Templates for the solution of large sparse linear systems consists of a collection of iterative methods together with a manual for algorithmic choices, instructions, and guidelines [6]. In contrast to the dense matrix case, there is no single iterative method that can solve any given sparse linear system in reasonable time and with reasonable memory requirements. In principle, many of the described methods solve a given system in at most iteration steps, where is the order of the linear system. Not only does this assume exact arithmetic, also one does not want to do as many as iteration steps for large systems. A major problem with iterative solution methods is that their speed of convergence depends on certain properties of the linear system, and these properties are often unknown in practice. This makes it difficult to make the right selection for a method. For many very large linear systems iterative systems are at this moment the only possible choice; there is hardly an alternative. Therefore, it is for engineers and scientists necessary to have a meaningful variety of iterative methods available, and this is what Templates provides. The potential user is guided into the world of iterative methods, and only a minimum of linear algebra knowledge is assumed. With the information provided in [6] the user can build high quality software at several levels of sophistication, depending on requirements and goals.

The starting point in Templates is the presentation of a selection of iteration methods in a very simple form, which closely follows the algorithmic descriptions in the original references. There is no monitoring of convergence behavior, and there are no precautions against numerical problems.

In this form the iteration schemes consist of one or two matrix vector products, a couple of inner-products and a few vector updates. The algorithms also allow for some form of preconditioning, which means that in each step a linear system has to be solved, where has to be specified by the user. All these operations are pretty well understood, and it is easy to see how the algorithms work. Although possible numerical complications are ignored, these basic iteration Templates may be used for simple computational experiments in order to get familiarized with the basic iteration methods and their mutual differences. The user is given extensive advice in the manual on how to expand the basic codes to reliable software for practical problems.

We believe that Templates in its present form is a good move in the process of educating and helping non-experts in understanding, and in the usage of iterative methods. Since the formulation is intentionally kept simple, it reveals the algorithmic essentials, and it shows where the major computational effort is going to be spent. The Matlab codes facilitate the user to play with the algorithms for relatively small but representative linear systems on workstations or even PC's, and help to develop some feeling for the potential usefulness of each of the methods. Together with the manual the chosen formulations provide the non-expert with a low-level entrance to the field. This approach has received approval recently in some new textbooks [40][39].

The formulations and the notations have been chosen so that differences and similarities between algorithms are easily recognized. In this respect they provide simple and easy to code algorithms, that can serve as a yardstick with which other implementations can be compared. Templates provides codes in Fortran, Matlab, and C. The codes are equivalent for each algorithm, and the only differences in outputs may be attributed to the respective computer arithmetic. This usage of Templates makes it easier to compare and to present results for iteration methods in publications. It is easier now to check results which have been obtained through these standardized codes.

As an example of how Templates is intended to function, we discuss the algorithm for the preconditioned Conjugate Gradient method for the iterative solution of the linear system , with a symmetric positive definite matrix:

In order to make the code useful for a specific problem class, the user has to provide a subroutine that computes , for a given vector . In Templates such code is not specified, but the user is given instructions on how to exploit sparsity of the given matrix [6], and by well-chosen examples it is shown how to code the matrix vector product for specific data storage formats. The following sparse matrix storage formats are discussed at some length:

1. Compressed Row Storage
2. Compressed Column Storage
3. Block Compressed Row Storage
4. Compressed Diagonal Storage
5. Jagged Diagonal Storage
6. Skyline Storage
For example, in Compressed Row Storage format subsequent nonzeros of the matrix rows are put in contiguous memory locations. In this case we create vectors: one for the floating point nonzero entries of (val), the other two with integers (col_ind, row_ptr). The row_ptr vector stores the locations in the val vector that start a row; the col_ind vector stores the column indexes of the elements in the val vector. The matrix vector product should then look like

Also suggestions are supplied for coding this on parallel computers.
For some algorithms also code has to be provided that generates , and this is given much attention.

Similarly, the user has to provide a subroutine that solves the equation for given right-hand side . In this respect is called the preconditioner and the choice for an appropriate preconditioner is left to the user. The main reason for this is that it is very difficult, if not impossible, to predict what a good preconditioner might be for a given problem of which little else is known besides the fact that is symmetric positive definite. In [6] an overview is given of possible choices and semi-code for the construction of some popular preconditioners is provided, along with instructions and examples on how to code them for certain data-storage formats. A popular preconditioner that is often used in connection with the Conjugate Gradient method is the Incomplete Cholesky preconditioner. If we write the matrix as ( is strictly lower diagonal, is diagonal), then a simplified form of the preconditioner can be written as , and in [6] formulas are given for the computation of the elements of the diagonal matrix . Since the factor is part of , the coding for the preconditioner is not too complicated.

The other operations (inner-products and vector updates) in the algorithm are coded, in the Fortran codes, as BLAS operations. On most modern computers these BLAS are standard available as optimized kernels, and the user needs not to be concerned with these vector operations.

So if the user wants to have transportable code then he only needs to focus on the matrix vector product and the preconditioner. However, there is more that the user might want to reflect on. In the standard codes only a simple stopping criterion is provided, but in [6] a lengthy discussion is provided on different possible termination strategies and their merits. Templates for these strategies have not yet been provided since the better of these strategies require estimates for such quantities as , or even , and we do not know about convenient algorithms for estimating these quantities (except possibly for the symmetric positive definite case). This is in contrast with dense systems, for which efficient estimators are available. As soon as a reliable estimating mechanisms for these norms become known then a template can be added.

The situation for other methods, for instance Bi-CG, may be slightly more complicated. Of course the standard template is still simple but now the user may be facing (near) break-down situations in actual computing. The solution for this problem is quite complicated, and the average user may not be expected to incorporate this easily and is referred to more professional software at this moment. A template for a cure of the break-downs would be very nice but is not planned as yet.

Also the updating for the residual vector and the approximate solution may suffer from problems due to the use of finite precision arithmetic. This may result into a vector that differs significantly from , whereas these vectors should be equal in exact arithmetic. Recently, good and reliable strategies have been suggested [34], and these are easily incorporated into the present Templates.

At the moment the linear systems templates for iterative methods form a good starting point for the user for the construction of useful software modules. It is anticipated that in future updates, templates will be added for algorithmic details, such as reliable vector updating, stopping criteria, and error monitoring. It is also conceivable that templates will be provided that help the user to analyze the convergence behavior, and to help get more insight in characteristics of the given problem. Most of these aspects are now covered in the Templates book, in the form of guidelines and examples. The Templates book [6] gives an abundance of references for more background and for all kinds of relevant details.

Next: Templates for Solution Up: Templates for Linear Algebra Previous: The BLAS as

Jack Dongarra
Wed Jun 21 02:35:11 EDT 1995