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:
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:
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. |
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
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 |
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.
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.
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.
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 )After zeroing out NONZ, DIAG and ENV, store the numerical values of the matrix and the right hand side. You can use the user routines as shown below.
call addrqt ( isub, jsub, value, invprm, diag, xenv, env, xnonz, nonz, nzsubs, n ) call addrhs ( invprm, isub, n, rhs, value )Factor and solve the system.
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 )Unpermute the solution stored in RHS.
call perm_rv ( n, rhs, perm )
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 )Zero out the matrix and right hand side storage in DIAG, LNZ, and RHS, then store numerical values. You can build up the values by calling the following routines repeatedly.
call addcom ( isub, jsub, value, invprm, diag, lnz, xlnz, nzsub, xnzsub, n ) call addrhs ( invprm, isub, n, rhs, value )Factor and solve the system.
call gs_factor ( n, xlnz, lnz, xnzsub, nzsub, diag, link, first ) call gs_solve ( n, xlnz, lnz, xnzsub, nzsub, diag, rhs )Unpermute the solution stored in RHS.
call perm_rv ( n, rhs, perm )
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 )Zero out the matrix and right hand side storage in DIAG, LNZ and RHS. Store numerical values, perhaps using the following routines repeatedly.
call addcom ( isub, jsub, value, invprm, diag, lnz, xlnz, nzsub, xnzsub ) call addrhs ( invprm, isub, n, rhs, value )Factor and solve.
call gs_factor ( n, xlnz, lnz, xnzsub, nzsub, diag, link, first ) call gs_solve ( n, xlnz, lnz, xnzsub, nzsub, diag, rhs )Unpermute the solution stored in RHS.
call perm_rv ( n, rhs, perm )
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 )Now store actual matrix values in DIAG, ENV and RHS, calling the following two routines as often as necessary. Be sure to first zero out the storage for DIAG, ENV and RHS. Only call ADDRCM for A(I,J) with I greater than or equal to J, since only the lower half of the matrix is stored.
call addrcm ( isub, jsub, value, invprm, diag, xenv, env, n ) call addrhs ( invprm, isub, n, rhs, value )Now factor and solve the system.
call es_factor ( n, xenv, env, diag ) call el_solve ( n, xenv, env, diag, rhs ) call eu_solve ( n, xenv, env, diag, rhs )Unpermute the solution stored in RHS.
call perm_rv ( n, rhs, perm )
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 )After zeroing out NONZ, DIAG and ENV, store the values of the matrix and the right hand side. You can use the user utility routines as shown below.
call addrqt ( isub, jsub, value, invprm, diag, xenv, env, xnonz, nonz, nzsubs, n ) call addrhs ( invprm, isub, n, rhs, value )Factor and solve the system.
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 )Unpermute the solution stored in RHS.
call perm_rv ( n, rhs, perm )
SPARSEPAK is available in a FORTRAN77 version and a FORTRAN90 version.
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.
You can go up one level to the FORTRAN77 source codes.