SPARSEPAK
(Obsolete) Waterloo Sparse Matrix Package

SPARSEPAK is a FORTRAN77 library which solves large sparse systems of linear equations.

SPARSEPAK is an old version of the Waterloo Sparse Matrix Package. SPARSEPAK 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.

SPARSPAK 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 )
      

Languages:

SPARSEPAK is available in a FORTRAN77 version and a FORTRAN90 version.

Related Data and Programs:

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

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

CSPARSE, a C library which contains direct sparse matrix operations.

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;

Examples and Tests:

• sparsepak_prb.f, a sample problem;
• sparsepak_prb.sh, commands to compile, link and run the sample problem;
• sparsepak_prb_output.txt, output from the sample problem;

List of Routines:

• ADDCOM adds values to a matrix stored in compressed storage scheme.
• ADDRCM adds values to a matrix stored in the RCM scheme.
• ADDRHS adds a quantity to a specific entry of the right hand side.
• ADDRQT adds values to a matrix stored in the implicit block storage scheme.
• ADJ_ENV_SIZE computes the envelope size for an adjacency structure.
• ADJ_PRINT prints the adjacency information stored in ADJ_ROW and ADJ.
• ADJ_SET sets up the adjacency information ADJ_ROW and ADJ.
• ADJ_SHOW displays a symbolic picture of a matrix.
• BLOCK_SHUFFLE renumbers the nodes of each block to reduce its envelope.
• DEGREE computes node degrees in a connected component, for the RCM method.
• EL_SOLVE solves a lower triangular system stored in the envelope format.
• ES_FACTOR factors a positive definite envelope matrix into L * L'.
• EU_SOLVE solves an upper triangular system stored in the envelope format.
• FNBENV finds the envelope of the diagonal blocks of a Cholesky factor.
• FNDSEP finds a small separator for a connected component in a graph.
• FNENV finds the envelope structure of a permuted matrix.
• FNLVLS generates a rooted level structure, as part of the RQT method.
• FNOFNZ finds columns of off-block-diagonal nonzeros.
• FNSPAN finds the span of a level subgraph subset, as part of the RQT method.
• FNTADJ determines the quotient tree adjacency structure for a graph.
• FNTENV determines the envelope index vector.
• FN1WD finds one-way dissectors of a connected component for the 1WD method.
• GENND finds a nested dissection ordering for a general graph.
• GENQMD implements the quotient minimum degree algorithm.
• GENRCM finds the reverse Cuthill-Mckee ordering for a general graph.
• GENRQT determines a partitioned ordering using refined quotient trees.
• GEN1WD partitions a general graph, for the 1WD method.
• GS_FACTOR: symmetric factorization for a general sparse system.
• GS_SOLVE solves a factored general sparse system.
• I4_SWAP switches two integer values.
• I4VEC_COPY copies one integer vector into another.
• I4VEC_INDICATOR sets an integer vector to the indicator vector.
• I4VEC_REVERSE reverses the elements of an integer vector.
• I4VEC_SORT_INSERT_A uses an ascending insertion sort on an integer vector.
• LEVEL_SET generates the connected level structure rooted at a given node.
• PERM_INVERSE produces the inverse of a given permutation.
• PERM_RV undoes the permutation of the right hand side.
• QMDMRG merges indistinguishable nodes for the QMD method.
• QMDQT performs the quotient graph transformation.
• QMDRCH determines the reachable set of a node, for the QMD method.
• QMDUPD updates the node degrees, for the QMD method.
• RCM renumbers a connected component by the reverse Cuthill McKee algorithm.
• RCM_SUB finds the reverse Cuthill McKee ordering for a given subgraph.
• REACH determines the reachable set of a node through a subset in a subgraph.
• ROOT_FIND finds pseudo-peripheral nodes.
• RQTREE finds a quotient tree ordering for a component, in the RQT method.
• SMB_FACTOR performs symbolic factorization on a permuted linear system.
• TIMESTAMP prints the current YMDHMS date as a time stamp.
• TS_FACTOR performs the symmetric factorization of a tree-partitioned system.
• TS_SOLVE solves a tree-partitioned factored system.

You can go up one level to the FORTRAN77 source codes.