Comments and Article Updates
 

Chapter 1: Introduction

Chapter 2: Discussion

2.1   Preconditioners
2.2   Data structure issues
2.3   High performance computing
2.4   Language

3.1   Aztec
3.2   BlockSolve95
3.3   BPKIT2
3.4   GPS: General Purpose Solver
3.5   IML++
3.6   Itpack 2C / ItpackV 2D
3.7   Laspack
3.8   ParPre
3.9   PCG
3.10 Petsc
3.11 Parallel Iterative Methods (PIM)
3.12 PSparselib
3.13 Sparse Approximate Inverse (SPAI) Preconditioner
3.14 SPlib
3.15 Templates

Chapter 4: Comparison Chart of Features

Chapter 5: Performance tests

5.1 Machines used
5.2 Results

Chapter 1 Introduction

Table of Contents

NHSE Review Home

Chapter 2 Discussion

2.1 Preconditioners

2.1.1 About incomplete factorisation preconditioners

BlockSolve95 is the only package catching breakdown of the ILU or IC factorisation. The manual outlines code that successively shifts the diagonal until the factorisation succeeds.

2.1.2 Preconditioners for sequential processing

2.1.3 Preconditioners for parallel iterative methods

Direct approximations of the inverse SPAI (section 3.13) is the only package that provides a direct approximation method to the inverse of the coefficient matrix. Since such an approximation is applied directly, using a matrix-vector product, it is trivially parallel. The SPAI preconditioner is in addition also generated fully in parallel.

Block Jacobi Each processor solves its own subsystem, and there is no communication between the processors. Since this strategy neglects the global/implicit properties of the linear system, only a limited improvement in the number of iterations can result. On the other hand, this type of method is very parallel.

All parallel preconditioner packages provide some form of Block Jacobi method.

Additive Schwarz As in the Block Jacobi method, each processor solves a local subsystem. However, the local system is now augmented to include bordering variables, which belong to other processors. A certain amount of communication is now necessary, and the convergence improvement can by much higher.

This method is available in Aztec (3.1), Petsc (3.10), ParPre (3.8), PSparselib (3.12).

Multicolour factorisations It is possible to perform global incomplete factorisation if the domain is not only split over the processors, but is also augmented with a multicolour structure. Under the reasonable assumption that each processor has variables of every colour, both the factorisation and the solution with a system thus ordered are parallel. The number of
synchronisation points is equal to the number of colours.

This is the method supplied in BlockSolve95 (3.2); it is also available in ParPre (3.8).
 

Block factorisation It is possible to implement block SSOR or ILU methods, with the subblocks corresponding to the local systems (with overlap, this gives the Multiplicative Schwarz method). Such factorisations are necessarily more sequential than Block Jacobi or Additive Schwarz methods, but they are also more accurate. With an appropriately chosen processor ordering (e.g., multicolour instead of sequential) the solution time is only a small multiple times that of Block Jacobi and Additive Schwarz methods.

Such block factorisations are available in Parpre (3.8); PSparselib (3.12) has the Multiplicative Schwarz method, under the name `multicolour SOR'.
 

Multilevel methods Petsc (3.10) and ParPre (3.8) are the only packages supplying variants of (algebraic) multigrid methods in parallel.

Schur complement methods ParPre (3.8) and PSparselib (3.12) contain Schur complement domain decomposition methods.

2.2 Data structure issues

2.2.1 The interface to user data

The packages reviewed here have taken the following list of approaches to this problem.
 
 

• Fully internal, not directly accessible, data structures (basically object oriented programming). This approach is taken in Petsc (3.10), ParPre (3.8),and BlockSolve (3.2).  The data structures here are only accessible through function calls. This means that the package has to supply routines for constructing the structures, for inspecting them, and for applying them.
• Prescribed, user supplied, data structures. This approach is taken in Aztec (3.1); it is available in PCG (3.9).Here the user has to construct the data structures, but the package supplies the product and solve routines.
• User supplied, arbitrary, data structures, with user supplied product and solve routines. This approach is available in PCG (3.9) and Petsc (3.10); the object-oriented package IML++ (3.5) can also be considered to use this approach.
• User defined data structures, not passed to the iterative method at all; product and solve operations are requested through a reverse communication interface. This approach is taken in PIM (3.11); it is also available in PCG (3.9).

2.2.2 Parallel data layout

• Each processor receives a contiguous series of rows of the matrix. This is the approach taken in Petsc (3.10); it is available in BlockSolve95 (3.2). Under this approach a processor can determine the ownership of any variable, by keeping a short list of first and last rows of each processor.
• Each processor receives an arbitrary set or rows. This is the approach taken in Aztec (3.1); it is available in BlockSolve95 (3.2). Under this scheme, each processor needs to maintain the mapping of its local rows to the global numbering. It is now no longer trivial to determine ownership of a variable.

2.3 High performance computing

Other performance problems stem from parallel aspects of the solution methods.

Here are then a few of the approaches taken to alleviate performance problems.
 

2.3.1 Quotient graph computation

However, sparse matrices from problems with multiple physical variables will have dense subblocks, for instance corresponding to the variables on any given node. It is possible to increase the efficiency of linear algebra operations by imposing on the matrix the block structure induced by these small dense subblocks. For instance, the scalar multiplication/division a_ika_kk^-1a_ki that appears in Gaussian elimination now becomes a matrix inversion and a matrix-matrix multiplication, in other words, BLAS3 kernels.

Identifying cliques in the matrix graph, and deriving a quotient graph by `factoring them out', is the approach taken in the BlockSolve package; section 3.2.

The PCG package (section 3.9) takes the opposite approach in its Regular Grid Stencil data format. To get dense subblocks, the degrees of freedom have to be numbered first in regular storage, but in PCG they are numbered last. This approach is more geared towards vector architectures.

2.3.2 Inner products

The conjugate gradients and bi-conjugate gradients algorithms have two interdependent inner products, and various people have suggested ways to reduce this to two independent ones, which can then combine their communication stages. See [3] and the references cited therein. This approach has been taken in the PIM package; section 3.11.

The second source of inner products is the orthogonalisation stage in the GMRES algorithm. Using Gram-Schmidt orthogonalisation, the inner products are independent and can be combined; however, this makes the resulting algorithm unstable. Modified Gram-Schmidt gives a more stable algorithm, but the inner products are interdependent, and have to be processed in sequence. A middle ground is to apply (unmodified) Gram-Schmidt twice.
 

2.4 Language

Table of Contents
 

NHSE Review Home

Chapter 3 The Packages

3.1 Aztec

3.1.1 Basic information

3.1.2 Contents

Preconditioners Block Jacobi with ILU on the subblocks; also Additive Schwarz, with an overlap limited to one.

Data structures Distributed variants of Compressed Row and Variable Block Row; see below.

Manual User's Guide, 43 pages

Example codes Yes
 

3.1.3 Parallelism and data layout

3.1.4 Other

3.1.5 Installation

•  I do not see any advantage over the CRS format.
• While conversion to and from CRS is straightforward, the arrays need to be longer by one position. This means that interfacing Aztec to a code using CRS entails deallocating and reallocating the arrays.
• The extra storage location in the real array is not used; the location in the integer array duplicates the scalar parameter giving the number of rows.

3.2 BlockSolve95

BlockSolve can be used with the Petsc package; section 3.10.

3.2.1 Basic information

3.2.2 Contents

Preconditioners diagonal scaling, block Jacobi with blocks corresponding to the cliques factored out of the quotient graph, incomplete LU and Cholesky.

Data structures Internal, as a C structure..
 
 
 
 

Example codes Yes

 3.2.3 Parallelism and data layout

1. Since the definition of the data structure is given in the manual, the user can explicitly construct it.
2. The user can supply a compressed row storage matrix, with column indices in the local numbering (section 2.2.2), to the routine BSeasy_A, which yields the matrix structure.

3.3 BPKIT2

3.3.1 Basic information

3.3.2 Contents

Preconditioners Block SSOR and ILU, possibly with block fill-in, and with various methods for explicit and implicit approximation of inverses of pivot blocks

Data structures Internal, as a C++ class.
 
 
 
 

Example codes Yes, in Fortran, C, and C++
 

3.3.3 Installation

The manual is not well written. Many parameters and auxiliary routines are under-documented.
 

3.4 GPS: General Purpose Solver

3.4.1 Basic information

3.5 IML++

An example implementation of such classes called Sparselib++ is available from the same authors. It contains uni-processor matrix classes based on compressed row and column and coordinate storage, a vector class, and Jacobi and ILU/IC preconditioners for the compressed row and column matrices. Additionally it contains tools for converting a Harwell-Boeing matrix or matrix file to these formats.
 

3.5.1 Basic information

3.5.2 Contents

Preconditioners n/a

Data structures n/a Manual Reference guide, 39 pages; also SparseLib++ Reference Guide, 20 pages.

Example codes no
 
 
 

3.6 Itpack 2C / ItpackV 2D

Itpack features adaptive determination of matrix bounds, to be used in accurate stopping tests, of for the Chebyshev semi-iteration.
 

3.6.1 Basic information

3.6.2 Contents

Preconditioners n/a; see above

Data structures Itpack 2C: Compressed Row (there are auxiliary routines to facilitate building the data structure); ItpackV 2D: ellpack storage.

Manual 22/14 pages

Example codes Yes
 
 

3.6.3 Installation

The PROGRAM statement in the test files had to be edited.

The test code was insufficiently documented for an easy 'same problem but larger' test.

Nspcg.  The nspcg routines come in files with undescriptive names such as nspcg1.f.

The nspcg5.f file needed replacement of the timer routine.
 

3.7 Laspack

There are data structures for vectors, general and square matrices, and preconditioners; a large number of accessing and manipulating these objects is available.

3.7.1 Basic information

3.7.2 Contents

Multigrid methods Conventional and Full multigrid, BPX preconditioning

Preconditioners Jacobi, SSOR, ILU(0)

Data structures Internal, created by passing elements by function call.
 
 
 

Example codes Yes.

3.8 ParPre

ParPre can be used independently of Petsc, but does require Petsc to be installed.

3.8.1 Basic information

3.8.2 Contents

Preconditioners Additive and multiplicative Schwarz, Generalised Block SSOR, Schur complement domain decomposition, Algebraic multilevel methods (including multicolour ILU and algebraic multigrid).

Data structures Internal.
 
 

Example codes Yes
 

3.8.3 Parallelism and data layout

For the methods that do not have all subdomains operate in parallel (e.g., multiplicative Schwarz as opposed to additive Schwarz), the user can specify the strategy that determines the sequential ordering of the domains. The choices are: natural, red-black, and multicolour.

3.9 PCG

3.9.1 Basic information

3.9.2 Contents

Preconditioners Richardson and Jacobi

Data structures Regular Grid Stencil storage supported; also data structure free mode by passing product and solve routines, or through reverse communication.

Manual Reference manual, 64 pages; also Getting Started manual, Examples manual, and guide to the XPCG Xwindows front end.

Example codes Yes.

3.9.3 Interface to user data structures

First of all, PCG has one supported data structure: the Regular Grid Stencil storage. This corresponds to the type of matrices arising from having the same finite difference or finite element stencil on each grid point of a Cartesian product grid. Such matrices are trivially stored in multi-dimensional arrays. After the user sets up the matrix array and the right hand side vector, PCG solves the system (this is called Top Level Usage of PCG).
 

Secondly, the user can pass matrix and preconditioner arrays, plus two integer arrays for each, to iterative method routines on the Iterative Method Level of PCG. Additionally, now the user has to supply routines to perform the matrix vector product and the preconditioner solve, plus their transposes for methods such as BiCG that need them.

Thirdly, PCG supports Reverse Communication Usage: no matrix or preconditioner structures whatsoever are passed to the iterative routines. Instead, when a product or solve operations needs to be performed, the iterative routine will save its internal state, and return control to the user program with a request to perform that operation.

1     continue
      call CG( .... IREQ .... IVA, IVQL ... FWK ...)

      if ( IREQ .eq. JAV ) then
c perform matrix-vector product to the vector FWK(IVQR)
c leaving the result in FWK(IVA)
      else if ( IREQ .eq. JQLV )
c apply left preconditioner to the vector FWK(IVA)
c leaving the result in FWK(IVQR)
      else .....
      end if

      goto 1

Control is also returned to the user for the convergence test, but inner products are still performed inside the iterative routine.

The Regular Grid Stencil storage scheme incorporates the possibility of having multiple physical variables per grid point; see section 2.3.1.
 

3.10 Petsc

Petsc can use the preconditioners of the BlockSolve package; section 3.2.

3.10.1 Basic information

3.10.2 Contents

Preconditioners Identity, Jacobi, Block Jacobi, Block Gauss-Seidel (only sequential), SOR and SSOR, IC and ILU (sequential only?), Additive Schwarz, full factorisation (sequential only), user supplied.

Data structures Internal, elements passed through function calls.
 
 

Example codes Yes

3.10.3 Parallelism and data layout

The data structures are completely hidden from the user, only accessible through function calls. Creation of the vector and matrix data structures is
completely general: any processor can specify elements for any other processor.

3.11 Parallel Iterative Methods (PIM)

3.11.1 Basic information

3.11.2 Contents

Preconditioners None

Data structures None

Manual User's guide; 81 pages

Example codes Yes; both sequential and parallel, and for dense and sparse data formats. The manual contains a good discussion of the example programs.
 

3.11.3 Interface to user data structures

The calling interface for these external routines is fairly simple, e.g.,

subroutine matvec{u,v,ipar}
double precision u(*),v(*)
integer ipar(*)

where the ipar array is the information array that is passed to the iterative method.

Unfortunately this puts a number of restriction on the user's choices. For instance, it implies that the matrix has to be in common blocks, and that the vectors need to be simple arrays; they can not be pointers to more elaborate structures.
 

3.12 PSparselib

3.12.1 Basic information

Author(s) Yousef Saad and Gen-Ching Lo

Latest version 2.15, May 1997 (manual is dated June 1996)

Status Being developed; future version may be for internal use only.
 

3.12.2 Contents

Preconditioners Additive and multiplicative Schwarz, Schur complement domain decomposition

Data structures Undocumented, the user is referred to tech reports and articles.

Manual Users manual, 14 pages; this does not document calling sequences or data structures.

Example codes Yes
 

3.12.3 Parallelism and data layout

3.12.4 Installation

3.13 Sparse Approximate Inverse (SPAI) Preconditioner

Note:  SPAI's page no longer exists. The developers report problems with the current version of SPAI and are working on a much-improved version.
 

3.13.1 Basic information

3.13.2 Approximate inverse

The approximate inverse class of preconditioners is different in that they compute explicitly an approximation to the inverse. Hence the application is an explicit matrix-vector product operation, and therefore trivially parallel.

The method in the SPAI code is based on ideas from [5]: the minimisation problem
 
 
 

is solved, with the sparsity pattern of M predetermined or adaptively determined. This minimisation problem turns out to reduce to independent subproblems for the rows or columns of M, and is therefore executable in parallel.

An other advantage of this method is that it is not subject to breakdown the way factorisation based methods are.

3.14 SPlib

3.14.1 Basic information

3.14.2 Contents

Preconditioners Identity, ILU(s), MILU(s,r), ILUT(s,t), SSOR(omega), TRID(s), ILU0, ECIMGS; where s is the number of levels of fill, r is the relaxation parameter [1], t is the drop tolerance.

Data structures Compressed Sparse Row

Manual 26 pages

Example codes Driver program that read a Harwell-Boeing matrix and solves a problem with it.
 

3.15 Templates

3.15.1 Basic information

3.15.2 Contents

Preconditioners Identity, Jacobi

Data structures User supplied: each iterative method is given in two versions.
 
 
 

Example codes Yes

Table of Contents
 

NHSE Review Home

Chapter 4 Comparison chart of features

Parallel Does the package run in parallel? All the parallel packages are based on MPI, other protocols are noted.

Iterative Does the package contain iterative methods? A few packages have preconditioners as their main focus, but suppply one or a few iterative methods for the user who doesn't have any yet.

Prec Does the package contain preconditioners?

Data How does the package interface to user data? See note 3 below.

Lang What is the implementation language of the package?

Inst Is the library instrumented, reporting on flops and timing?
 
 
 
 
 
  

Footnotes
¹ Not the main focus of this package.

² The library abstracts away from data structure implementation aspects; parallelism is possible, but is the user's responsibility.

³ For the explanation of terms `internal', `prescribed', `supplied', and `free', see section 2.2.1.
^3a converted from compressed row format.
^3b converted from Harwell-Boeing format.

⁴ Identical to Petsc format.

⁵ Nothing beyond Jacobi.

⁶ The external product and solve routines are presumed to find the matrix in a common block.

⁷ Can not be chosen independently of the iterative method: the user picks a combination.

⁸ Fortran interface provided.

⁹ Preconditioners provided in an example C++ matrix class library, SparseLib++.
 
 

 Table of Contents

 NHSE Review Home
 

Chapter 5 Performance tests

• Single processor and parallel where possible.
• Using all data structures supplied or accepted by the package.
• Comparing various iterative methods, in particular where they have different parallel behaviour, such as the Chebyshev method versus the Conjugate Gradient method.
• Assessing the efficacy of different preconditioners, measuring separately
• and combined:

• Solving different systems of the same size and structure is not of much use, as this only changes the number of iterations performed; one could note how many iterations are necessary to offset the initial setup cost.

5.1 Machines used

nala Sun Ultra-Sparc 2200 with Solaris 5.5.1. 200MHz, 16K L1, 1Mb L2. Compilers: f77 -O5 and gcc -O2.

cetus lab Sun Ultra-Sparc 1, 143 Mhz, connected by 10Mbps Ethernet.
 
 
 
 
  

5.2 Results

5.2.1 Aztec

• CG instead of CGS,
• Block Jacobi preconditioner,
• 5-point instead of 7-point matrix.

Aztec's built in timing and flop count does not support the ILU(0) preconditioner, so we added that. The flop count is approximate, but does not overestimate by more than a few percent. We omit the N = 400 tests because they were too short for the timer to be reliable.

From table 2 we see for small problem sizes the communication overhead dominates; for larger problems the performance seems to level off at 13 Mfl per processors, about 10 percent of peak performance. Performance of an ILU(0)- preconditioned method (table 3) is slightly lower. The explanation for this is not immediately clear. Note that, since we used a regular grid problem, it is
not due to indirect addressing overhead.
 
 
 
 

5.2.2 BlockSolve95

From table 4 we see that the performance of BlockSolve95 is less than of other packages reported here. This is probably due to the more general data format and the resultant indirect addressing overhead. Results in table 5 show that by inode/clique identification BlockSolve95 can achieve performance comparable to regular grid problems in other packages.

Larger problems than those reported led to segmentation faults, probably because of insufficient memory. Occasionally, but not always, BlockSolve aborts with an 'Out of memory' error.
 
 

5.2.3 Itpack

Certain Itpack files are provided only in single precision. We took the single precision files and compiled them with f77 -r8 -i4, which makes the REALs 8 bytes and INTEGERs 4. It is not clear why diagonal storage will only give good performance on larger problems.
 
 

5.2.4 Petsc

We wrote our own main program to generate the five-point Laplacian matrix. The method in table 7 is an unpreconditioned CG algorithm.

We tested the efficacy of ILU by specifying

PCSetType(pc,PCSOR);
PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP);

5.2.5 PSparsLib

Larger problem sizes ran into what looks like a memory-overwrite. Attempts to allocate more storage failed.

5.3 Discussion

• A sequential code should be able to attain 10--20% of the peak speed of the machine. This value was attained by most packages, using a variety of storage formats.
• Parallel codes have significant overhead; for large enough problems this is amortized to where the per-processor performance is about half of that of the sequential code
• Inode/clique identification can make a large difference in systems that have multiple variables per node.

References

[2] Richard Barrett, Michael Berry, Tony F. Chan, James Demmel, June Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo, Charles Romine, and Henk van der Vorst. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. SIAM, Philadelphia PA, 1994.

[3] E.F. D'Azevedo, V.L. Eijkhout, and C.H. Romine. Lapack working note 56: Reducing communication costs in the conjugate gradient algorithm on distributed memory multiprocessor. Technical Report CS-93-185, Computer Science Department, University of Tennessee, Knoxville, 1993. to appear.

[4] Jack Dongarra and Henk van der Vorst. Performance of various computers using sparse linear equations software in a fortran environment. Supercomputer, 1992.

[5] L. Yu. Kolotilina and A. Yu. Yeremin. On a family of two-level preconditionings of the incomplete block factorization type. Sov. J. Numer. Anal. Math. Modelling, pages 293--320, 1986.
 
 
 

Table of Contents
 NHSE Review Home

Copyright © 1998