Skyline Solution of a Sparse System of Linear Equations

The finite difference equations along with external and internal boundary conditions and storage area equations result in a system of linear equations which must be solved for each time step:

For a single channel without a storage area, the coefficient matrix has a band width of five and can be solved by one of many banded matrix solvers.

For network problems, sparse terms destroy the banded structure. The sparse terms enter and leave at the boundary equations and at the storage areas. The figure below shows a simple system with four reaches and a storage area off of reach 2.

The corresponding coefficient matrix is shown in the figure below. The elements are banded for the reaches but sparse elements appear at the reach boundaries and at the storage area. This small system is a trivial problem to solve, but systems with hundreds of cross sections and tens of reaches pose a major numerical problem because of the sparse terms. Even the largest computers cannot store the coefficient matrix for a moderately sized problem, furthermore, the computer time required to solve such a large matrix using Gaussian elimination would be very large. Because most of the elements are zero, a majority of computer time would be wasted.

Three practical solution schemes have been used to solve the sparse system of linear equations: Barkau (1985) used a front solver scheme to eliminate terms to the left of the diagonal and pointers to identify sparse columns to the right of the diagonal. Cunge et al. (1980) and Shaffranekk (1981) used recursive schemes to significantly reduce the size of the sparse coefficient matrix. Tucci (1978) and Chen and Simons (1979) used the skyline storage scheme (Bathe and Wilson, 1976) to store the coefficient matrix. The goal of these schemes is to more effectively store the coefficient matrix. The front solver and skyline methods identify and store only the significant elements. The recursive schemes are more elegant, significantly reducing the number of linear equations. All use Gaussian elimination to solve the simultaneous equations.

A front solver performs the reduction pass of Gauss elimination before equations are entered into a coefficient matrix. Hence, the coefficient matrix is upper triangular. To further reduce storage, Barkau (1985) proposed indexing sparse columns to the right of the band, thus, only the band and the sparse terms were stored. Since row and column operations were minimized, the procedure should be as fast if not faster than any of the other procedures. But, the procedure could not be readily adapted to a wide variety of problems because of the way that the sparse terms were indexed. Hence, the program needed to be re-dimensioned and recompiled for each new problem.

The recursive schemes are ingenious. Cunge credits the initial application to Friazinov (1970). Cunge's scheme and Schaffranek's scheme are similar in approach but differ greatly in efficiency. Through recursive upward and downward passes, each single routing reach is transformed into two transfer equations which relate the stages and flows at the upstream and downstream boundaries. Cunge substitutes the transfer equations in which M is the number of junctions. Schraffranek combines the transfer equations with the boundary equations, resulting in a system of 4N equations in which N is the number of individual reaches. The coefficient matrix is sparse, but the degree is much less than the original system.

By using recursion, the algorithms minimize row and column operations. The key to the algorithm's speed is the solution of a reduced linear equation set. For smaller problems Gaussian elimination on the full matrix would suffice. For larger problems, some type of sparse matrix solver must be used, primarily to reduce the number of elementary operations. Consider, for example, a system of 50 reaches. Schaffranek's matrix would be 200 X 200 and Cunge's matrix would be 50 X 50, 2.7 million and 42,000 operations respectively (the number of operations is approximately 1/3 n³ where n is the number of rows).

Another disadvantage of the recursive scheme is adaptability. Lateral weirs which discharge into storage areas or which discharge into other reaches disrupt the recursion algorithm. These weirs may span a short distance or they may span an entire reach. The recursion algorithm, as presented in the above references, will not work for this problem. The algorithm can be adapted, but no documentation has yet been published.

Skyline is the name of a storage algorithm for a sparse matrix. In any sparse matrix, the non-zero elements from the linear system and from the Gaussian elimination procedure are to the left of the diagonal and in a column above the diagonal. This structure is shown in the figure below. Skyline stores these inverted "L shaped" structures in a vector, keeping the total storage at a minimum. Elements in skyline storage are accessed by row and column numbers. Elements outside the "L" are returned as zero, hence the skyline matrix functions exactly as the original matrix. Skyline storage can be adapted to any problem.

The efficiency of Gaussian elimination depends on the number of pointers into skyline storage. Tucci (1978) and Chen and Simons (1979) used the original algorithm as proposed by Bathe and Wilson (1976). This algorithm used only two pointers, the left limit and the upper limit of the "L", thus, a large number of unnecessary elementary operations are performed on zero elements and in searching for rows to reduce. Their solution was acceptable for small problems, but clearly deficient for large problems. Using additional pointers reduces the number of superfluous calculations. If the pointers identify all the sparse columns to the right of the diagonal, then the number of operations is minimized and the performance is similar to the front solver algorithm.

Skyline Solution Algorithm

The skyline storage algorithm was chosen to store the coefficient matrix. The Gauss elimination algorithm of Bathe and Wilson was abandoned because of its poor efficiency. Instead a modified algorithm with seven pointers was developed. The pointers are:

1. IDIA(IROW) - index of the diagonal element in row IROW in skyline storage.
2. ILEFT(IROW) - number of columns to the left of the diagonal.
3. IHIGH(IROW) - number of rows above the diagonal.
4. IRIGHT(IROW) - number of columns in the principal band to the right of the diagonal.
5. ISPCOL(J,IROW) - pointer to sparse columns to the right of the principal band.
6. IZSA(IS) - the row number of storage area IS.
7. IROWZ(N) - the row number of the continuity equation for segment N.

The pointers eliminate the meaningless operations on zero elements. This code is specifically designed for flood routing through a full network.