sparsepak

sparsepak, a Fortran77 code which solves large sparse systems of linear equations.

This is an old version of the Waterloo Sparse Matrix Package. It can carry out direct solution of large sparse linear systems. Only positive definite matrices should be used with this program. Three different storage methods, and five solution algorithms are supported.

The methods available are:

• 1WD: One way dissection, partitioned tree storage.
• ND: Nested dissection, compressed storage.
• QMD: Quotient minimum degree, compressed storage.
• RCM: Reverse Cuthill-McKee with envelope storage.
• RQT: Refined quotient tree, partitioned tree storage.

This version is not the most recent one. The most recent version requires a license from the University of Waterloo, and includes support for large sparse linear least squares problems and other features.

The code offers five different methods for the direct solution of a symmetric positive definite sparse linear system. Subroutines are provided which try to reorder the matrix before factorization so that as little storage as possible is required.

The five algorithms all follow the same steps:

1. Determine the structure of the matrix.
2. Reorder rows and columns for efficiency.
3. Store numerical values in matrix and right hand side.
4. Factor the matrix.
5. Solve the system (may be repeated for multiple right hand sides).
6. Unpermute the solution.

Methods available are:

Symbol	Description
1WD	One way dissection, partitioned tree storage.
ND	Nested dissection, compressed storage.
QMD	Quotient minimum degree, compressed storage.
RCM	Reverse Cuthill-McKee with envelope storage.
RQT	Refined quotient tree, partitioned tree storage.

Data Structures

The following is a list of the vectors required for the various methods, including their type, size, and whether they are set by the user or not. Unfortunately, for most of the arrays, although their values are determined by the package, their minimum size is not known beforehand, and is often hard to estimate. You must examine the reference to get a better idea of how large to dimension arrays like XBLK and NZSUB, for example.

          Values                             Used by:
Name      Set by    Size      Type      1WD  RCM  RQT  ND   QMD

IADJ      User      varies    I         X    X    X    X    X
BNUM                N         I         X         X
DEG                 N         I                             X
DIAG      User      N         R         X    X    X    X    X
ENV       User      varies    R         X    X    X
FATHER              NBLKS     I         X         X
FIRST               N         I         X         X    X    X
INVPRM              N         I         X    X    X    X    X
IXLNZ               N+1
LINK                N         I                        X    X
LNZ                           R                        X    X
LS                  N         I         X         X    X
MARKER              N         I         X              X    X
MRGLNK              N         I                        X    X
NBRHD                         I                             X
NODLVL              N         I                   X
NONZ                          R         X         X
NZSUB                         I                        X    X
NZSUBS                        I         X         X
OVRLP                                                       X
PERM                N         I         X    X    X    X    X
QLINK               N         I                             X
QSIZE               N         I                             X
RCHLNK              N         I                        X    X
RCHSET              N         I         X                   X
RHS       User      N         R         X    X    X    X    X
SEP                           I                        X
SUBG                N         I         X         X
TEMP                N         R         X         X    X    X
XADJ      User      N+1       I         X    X    X    X    X
XBLK                NBLKS+1   I         X         X
XENV      User      N+1       I         X    X    X
XLNZ                N         I                        X    X
XLS                 N+1       I         X         X    X
XNONZ               N+1       I         X         X
XNZSUB                        I                        X    X

Adjacency Information

You are responsible for defining the adjacency or nonzero structure of the original matrix. This information is required by the routines to efficiently reorder the system.

Suppose our linear system is a 6 by 6 matrix, where the entries marked 'X' are nonzero:

        XX000X
        XXXX00
        0XX0X0
        0X0X00
        00X0XX
        X000XX
      

The adjacency information records, for each row I, which columns (besides column I) contain a nonzero entry, and, for each column I, which rows besides row I are nonzero.

In the example, the nonzero list for row I=1 is 2 and 6. This is because A(1,2) and A(1,6) are nonzero. The whole list is

Row	Nonzero columns
1	2, 6
2	1, 3, 4
3	2, 5
4	2
5	3, 6
6	1, 5

Note that in the list of nonzeroes for row I, we never include I itself. It is assumed that I is nonzero, or will become so during the factorization.

The adjacency information consists of two vectors, IADJ and XADJ. The vector IADJ contains the concatenated list of nonzeroes. The auxilliary vector XADJ holds the location of the start of each row list in IADJ. Thus, XADJ(1) should always be 1. XADJ(2) should be the location in IADJ where the nonzero information for row 2 begins. Finally, XADJ(N+1) points to the first unused entry in IADJ. For the example,

Index	1	2	3	4	5	6	7	8	9	10	11	12	13
IADJ	2	6	1	3	4	2	5	2	3	6	1	3
XADJ	1	3	6	8	9	11	13

You can see that to retrieve information about row 3, you would first check XADJ(3) and XADJ(3+1), which have the values 6 and 8. This tells us that the list for row 3 begins at IADJ(6) and stops at IADJ(8-1). So the nonzeroes in row 3 occur in columns IADJ(6)=2 and IADJ(7)=5. This seems a lot of work, but this construction can be automated, and the storage savings can be enormous.

Envelope Information

Assume, as usual, that we are working with a symmetric matrix. The lower row bandwidth of row I, denoted by BETA(I), is defined in the obvious way to mean the distance between the diagonal entry in column I and the first nonzero entry in row I.

The envelope of a matrix is then the collection of all subdiagonal matrix positions (I,J) which lie within the lower row bandwidth of their given row. This concept is sometimes also called profile or skyline storage, particularly when columns are considered rather than rows.

To store a matrix in envelope format, we require two objects, an array ENV to hold the actual numbers, and a pointer array XENV to tell us where each row starts and ends.

Let us return to our earlier example:

        XX000X
        XXXX00
        0XX0X0
        0X0X00
        00X0XX
        X000XX
      

We can analyze this matrix as follows:

Row	Bandwidth
1	0
2	1
3	1
4	2
5	2
6	5

This tells us that there are a total of 11 entries in the lower envelope. XENV(I) tells us where the first entry of row I is stored, and XENV(I+1)-1 tells us where the last one is. For our example, then

Index	XENV
1	1
2	1
3	2
4	3
5	5
6	7
7	12

Note that the "extra" entry in XENV allows us to treat the last row in the same way as the other ones. Now the ENV array will require 11 entries, and they can be stored and accessed using the XENV array.

Setting Matrix Values

Once you have defined the adjacency structure, you have set up empty slots for your matrix. Into each position which you have set up, you can store a numerical value A(I,J). You can use the user utility routines ADDRCM, ADDRQT and ADDCOM to store values in A for you.

Alternatively, if you want to do this on your own, you must be familiar with the storage method used. For example, in the envelope storage scheme you must know, given I and J, in what entry of DIAG and ENV you should store the value. Interested readers are referred to the reference for more information on this.

Example of 1WD Method

Create the adjacency structure XADJ, IADJ, and reorder the variables and equations.

        call gen1wd ( n, xadj, iadj, nblks, xblk, perm, xls, ls )

        call perm_inverse ( n, perm, invprm )

        call fnbenv ( xadj, iadj, perm, invprm, nblks, xblk, xenv, ienv,
          marker, rchset, n )
      

        call addrqt ( isub, jsub, value, invprm, diag, xenv, env, xnonz, nonz,
          nzsubs, n )

        call addrhs ( invprm, isub, n, rhs, value )
      

        call ts_factor ( nblks, xblk, father, diag, xenv, env, xnonz, nonz,
          nzsubs, temp, first, n )

        call ts_solve ( nblks, xblk, diag, xenv, env, xnonz, nonz, nzsubs,
          rhs, temp, n )
      

        call perm_rv ( n, rhs, perm )
      

Example of ND Method

Create the adjacency structure XADJ, IADJ. Then reorder the variables and equations.

        call gennd ( n, xadj, iadj, perm, xls, ls )

        call perm_inverse ( n, perm, invprm )

        call smb_factor ( n, xadj, iadj, perm, invprm, xlnz, nofnz, xnzsub,
          nzsub, nofsub, rchlnk, mrglnk )
      

        call addcom ( isub, jsub, value, invprm, diag, lnz, xlnz, nzsub,
          xnzsub, n )

        call addrhs ( invprm, isub, n, rhs, value )
      

        call gs_factor ( n, xlnz, lnz, xnzsub, nzsub, diag, link, first )

        call gs_solve ( n, xlnz, lnz, xnzsub, nzsub, diag, rhs )
      

        call perm_rv ( n, rhs, perm )
      

Example of QMD Method

Create the adjacency structure in XADJ, IADJ. Reorder the variables and equations.

        call genqmd ( n, xadj, adjnc2, perm, invprm, marker, rchset,
          nbrhd, qsize, qlink, nofsub )

        call perm_inverse ( n, perm, invprm )

        call smb_factor ( n, xadj, iadj, perm, invprm, xlnz, nofnz, xnzsub,
          nzsub, nofsub, rchlnk, mrglnk )
      

        call addcom ( isub, jsub, value, invprm, diag, lnz, xlnz, nzsub, xnzsub )

        call addrhs ( invprm, isub, n, rhs, value )
      

        call gs_factor ( n, xlnz, lnz, xnzsub, nzsub, diag, link, first )

        call gs_solve ( n, xlnz, lnz, xnzsub, nzsub, diag, rhs )
      

        call perm_rv ( n, rhs, perm )
      

Example of RCM Method

Create the adjacency structure XADJ, IADJ. Then call the following routines to reorder the system.

        call genrcm ( n, xadj, iadj, perm, xls )

        call perm_inverse ( n, perm, invprm )

        call fnenv ( n, xadj, iadj, perm, invprm, xenv, ienv, bandw )
      

        call addrcm ( isub, jsub, value, invprm, diag, xenv, env, n )

        call addrhs ( invprm, isub, n, rhs, value )
      

        call es_factor ( n, xenv, env, diag )

        call el_solve ( n, xenv, env, diag, rhs )

        call eu_solve ( n, xenv, env, diag, rhs )
      

        call perm_rv ( n, rhs, perm )
      

Example of RQT Method

Create the adjacency structure XADJ, IADJ. Reorder the variables and equations.

        call genrqt ( n, xadj, iadj, nblks, xblk, perm, xls, ls, nodlvl )

        call block_shuffle ( xadj, iadj, perm, nblks, xblk, xls, n )

        call perm_inverse ( n, perm, invprm )

        call fntadj ( xadj, iadj, perm, nblks, xblk, father, n )

        call fntenv ( xadj, iadj, perm, invprm, nblks, xblk, xenv, ienv, n )

        call fnofnz ( xadj, iadj, perm, invprm, nblks, xblk, xnonz, nzsubs,
          nofnz, n )
      

        call addrqt ( isub, jsub, value, invprm, diag, xenv, env, xnonz, nonz,
          nzsubs, n )

        call addrhs ( invprm, isub, n, rhs, value )
      

        call ts_factor ( nblks, xblk, father, diag, xenv, env, xnonz, nonz,
          nzsubs, temp, first, n )

        call ts_solve ( nblks, xblk, diag, xenv, env, xnonz, nonz, nzsubs,
          rhs, temp, n )
      

        call perm_rv ( n, rhs, perm )
      

Licensing:

The computer code and data files described and made available on this web page are distributed under the MIT license.

Languages:

sparsepak is available in a Fortran77 version and a Fortran90 version.

Related Data and Programs:

sparsepak_test

ccs, a data directory which contains examples of the Compressed Column Storage (CCS) sparse matrix file format;

crs, a data directory which contains examples of the Compressed Row Storage (CRS) sparse matrix file format;

DLAP, a Fortran90 library which solves sparse linear systems.

HB_IO, a Fortran90 library which reads and writes sparse linear systems stored in the Harwell-Boeing Sparse Matrix format.

HB_TO_ST, a Fortran77 program which converts the sparse matrix information stored in a Harwell-Boeing file into a sparse triplet file.

MGMRES, a Fortran90 library which applies the restarted GMRES algorithm to solve a sparse linear system.

MM_IO, a Fortran90 library which reads and writes sparse linear systems stored in the Matrix Market format.

RCM, a Fortran90 library which applies reverse Cuthill-McKee reordering.

SPARSEKIT, a Fortran90 library which carries out various operations on sparse matrices, by Youcef Saad.

UMFPACK, a Fortran77 library which solves unsymmetric sparse linear systems, by Timothy Davis, Iain Duff.

Reference:

1. Alan George, Joseph Liu,
   Computer Solution of Large Sparse Positive Definite Matrices,
   Prentice Hall, 1981,
   ISBN: 0131652745,
   LC: QA188.G46.
2. Norman Gibbs, William Poole, Paul Stockmeyer,
   An Algorithm for Reducing the Bandwidth and Profile of a Sparse Matrix,
   SIAM Journal on Numerical Analysis,
   Volume 13, Number 2, April 1976, pages 236-250.
3. Norman Gibbs,
   Algorithm 509: A Hybrid Profile Reduction Algorithm,
   ACM Transactions on Mathematical Software,
   Volume 2, Number 4, December 1976, pages 378-387.
4. Donald Kreher, Douglas Simpson,
   Combinatorial Algorithms,
   CRC Press, 1998,
   ISBN: 0-8493-3988-X,
   LC: QA164.K73.

Source Code:

• sparsepak.f, the source code;
• sparsepak.sh, commands to compile the source code;