Introduction

Using the Unit

Portability

Updates for November 1999/September 2002/August 2005/January 2007

The Author

Download the unit sparsolv.zip (320k)

Includes source, test programs, and all the instructions on this page.

Unit SparSolv is designed to solve sparse linear systems. These are equation systems
of the form:

A.X = R

where A is the LHS matrix of coefficients size NxN, R is the RHS vector of constants size N, and X is a vector of variables size N.

It is designed to save time and space where A is fairly large but contains few non-zero elements. A typical example would be the need to solve 1000 simultaneous linear equations involving 1000 variables, where each equation used, on average, only 6 variables. To solve this system using conventional, non-sparse, methods, one would need to allocate a 1000x1000 matrix using 8 MB of RAM (it is advisable to use double precision for larger linear systems). Even if the RAM was available, it would take a long time to solve such a large matrix.

Sparse methods save space and time by only storing information about non-zero elements. SparSolv maintains an individual linked list for each equation. Each node in such a list contains an element value (Double: 8 bytes), a column number (Integer: 4 bytes) and a pointer (4 bytes) to the next node. Each node uses 16 bytes, which is twice the memory used to store one number in a conventional matrix. So, pursuing the example above, the 6000 (6 x 1000) non-zero elements would require only 96 kB of memory. For SparSolv, the RHS and various work vectors raise this to 140 kB.

SparSolv solves the linear system by reducing the LHS matrix to triangular form, using Gaussian elimination. It chooses the order of pivot rows and columns in such a way as to minimize the creation of new non-zero LHS elements. Nevertheless, such creation is unavoidable, and so more memory is continually required as the solution proceeds. The program TEST1 (see below), which solves the 1000x1000 0.6% density matrix mentioned above eventually required 1.2 MB of RAM, nearly 9 times the initial requirement - although much less than would be needed by non-sparse methods.

The speed and memory requirements of the algorithm depend not only on the size and density of the LHS, but also on its structure: the way the non-zero elements are arranged. Real-world problems often have a structure which the algorithm is able to exploit. As in some other computer problems (eg, sorting), no one algorithm will be best for all problems. I recommend SparSolv for systems of up to 12000 equations with up to 200000 non-zeros and up to 5% density. Your mileage will vary !

The core code of SparSolv has been used by me and colleagues for a number of years to solve systems of non-linear equations for economic models. These may be represented as a matrix equation F(n,x) = 0, where n is a vector of N endogenous (model-deterimined) variables, x is a vector of exogenous (user-determined) variables and there are N equations. By drawing on contemporary economic data, we usually have one set of values for n and x which satisfy F. We wish to find new values for n corresponding to some other vector of exogenous (policy) variables y. We apply the Newton-Raphson recurrence:

J(n,y).dn - = - F(n,y).

Here J is the Jacobian of F with respect to n, and dn is the vector of changes to n which we hope will cause the RHS to go to zero. So J is our LHS, -F the RHS and we seek to find dn using SparSolv.

To get the best from SparSolv, you have to understand a little linear algebra. See:

- Introductory:

Numerical Recipes (Press, Teukolsky, Vetterling and Flannery) ISBN 0 521 43108 5 - Advanced:

Matrix Computations (Golub and Van Loan) ISBN 0 8018 3772 3 - Sparse:

Direct Methods for Sparse Matrices (Duff, Erisman and Reid) ISBN 0 521 43108 5

If you want to solve really big systems, you need professional quality routines such as can be found on NetLib (though see update below). But these are available in Fortran or C rather than Pascal. I have had good results creating a Fortran DLL using the Harwell routines (authors same as last book above) which was called by my Delphi 2.0 program. This route is really only practical if you are using 32-bit Pascal such as Delphi 2.0 or later.

The simplest way to learn how to use the unit may be to study and run the two example programs:

- Test1.dpr: Example of using routines with N = 500.
- Test2.dpr: Example of using routines with N = 4, showing some of the errors that could occur.
- TestMTX.dpr: Windows program which solves matrices stored on disk in MTX format.

The method of storing sparse matrices on disk, or in the non-SparSolv parts of your program, is up to you. The supplied SparSolv routines create and populate sparse matrices in a special structure, which you do not need to understand (but if you want to study how it is done, the source is there). The TestMTX example shows ONE way of storing data on disk, but you could devise another.

If run through Delphi 1, Test1 and Test2 will appear as WinCrt [console] text windows. For Delphi 2 and above, they run as console apps. Non-Delphi users should rename the files to Test1.pas and Test2.pas, and compile from the DOS command line.

Unit SparSolv makes available 5 Boolean Functions and 3 procedures, listed below. Each function returns True if operation was successfully completed - otherwise False. If False is returned, call the procedure GetErrorMessage to find the reason. It IS important that you check the function return values.

The first step is to call InitStruc to initialize the sparse matrix storage structure. Then, call AddLHS and AddRHS a number of times to set the values for the LHS and RHS. Then call Solve1 to solve the system. Then GetAnswer once for each variable. Finally call ReleaseStruc to free the memory that has been used.

- function InitStruc(NumEq : Word) : Boolean;

Creates and initializes sparse matrix structure - call this first. NumEq is number of equations/variables. - function AddLHS(ThisEqu, ThisVar : Word; ThisVal : Double) : Boolean;

Add an element to sparse matrix for equation ThisEqu and variable ThisVar; if such an entry already exists, ThisVal will be added to existing value. - function AddRHS(ThisEqu : Word; ThisVal : Double) : Boolean;

Set RHS for equation ThisEqu; if RHS has already been set, ThisVal will be added to existing value. YOU MUST SET THE RHS FOR EACH EQUATION, even if you set it to zero. - function Solve1 : Boolean;

Calculate solutions; sparse matrix is destroyed. - function GetAnswer(ThisVar : Word; var ThisVal : Double) : Boolean;

Read solution for variable ThisVar - probably called for each variable in turn. - procedure ReleaseStruc;

Releases remaining memory used by sparse matrix structure - call this last. - procedure GetErrorMsg(var S : String; var N1, N2, N3 : Word);

N1: error number; S: Error Description; N2, N3 : other possibly useful numbers. - procedure ShowMat;

Displays small sparse matrix.

Unit SparSolv also makes available 3 LongInt variables to assist in monitoring and controlling memory use:

- read SparMemUsed to find no of bytes of heap currently used by routines.
- read MaxMemUsed to find max no of bytes used so far.
- set MaxMemToUse: max no of bytes routines are allowed to use: -1 means no limit: default is 4MB (16-Bit) or 2GB (32-Bit)

If any of the functions listed above returns false, use GetErrorMessage to find what the problem was. Each error number is called from a unique line of code, so you can search to find where the error was caused. The errors will fall into four categories:

1: You forgot to call the routines in the right order:

Initstruc, AddLHS and AddRHS, Solve1, GetAnswer, ReleaseStruc.

2: You passed out-of-range equation and variable numbers to AddLHS, AddRHS or GetAnswer.

3: You ran out of memory.

4: Your matrix was singular. This means that your equation system does not uniquely
determine all the variables. A matrix is singular if:

(a) any row or any column is empty or consists only of zeros.

(b) any row can be expressed as a linear combination of some other rows. In that case
Gaussian elimination will lead to an empty row.

(c) Like (b) but for columns.

The supplied program Test2 demonstrates various types of singular matrix error.

For sparse matrices, we distinguish between structural singularity (row or col has no elements) and numerical singularity (vital elements are zero).

In its checking phase Solve1 may fail with 'Empty Row', 'Row without Variables' or 'Empty Col' messages. These indicate that rows/cols had no elements at all. While scaling the LHS, Solve1 may report 'All- Zero Row' or 'All-Zero Column'. This means that rows/cols had elements, but they were all set to 0.0. Next Solve1 tries to save time by identifying variables which occur only once. If there are 2 or more of these in the same equation, it fails with the 'Two Singles' message.

Finally, during the main Gaussian elimination phase, Solve1 again fails if the pivot row has no elements or only zero elements, with messages 'Structurally Singular' or 'Numerically Singular' respectively.

Very many singular matrix errors arise from either poor input data or a poorly conceived equation system. If you suspect that SparSolv is wrongly reporting a matrix to be singular, please solve the matrix using a different routine before complaining to me!

- It only allows one Right Hand Side (RHS) column. Furthermore, unlike LU routines, you cannot solve once, then backsolve several times for different RHS. This also rules out iterative improvement.
- You cannot save and reuse the pivot sequence - this feature could be added.
- It only allows one sparse structure at a time: You must call ReleaseStruc before calling InitStruc a second time.
- It will not solve singular matrices, even if
**some**(but not all) variables are uniquely determined !

As supplied the code will run on any 32-bit Delphi up to and including BDS2006 and Turbo Professional. To port to other dialects, you may need to modify various defines {$IFDEF xxxx} near the top of some files. The code worked with Free Pascal Compiler version 2.0.0 (2005). Now mainly of historical interest, SparSolv has been tested with the following 16-bit Pascal versions:

- 16-bit real mode Turbo Pascal 7.0 {$VER70}

tpc -b test1.dpr - 16-bit protected mode Borland Pascal 7.0 {$VER70} {$DPMI}

bpc -cp -b test1.dpr - 16-bit protected mode Borland Pascal for Windows {$VER70} {$WINDOWS}

bpc -cw -b test1.dpr - 16-bit Windows Delphi 1.0 {$VER80} {$WINDOWS}

dcc -b test1.dpr

An old-fashioned programming style has been used - no classes, objects, or modern syntactical conveniences like BREAK or CONTINUE. The reason is that when the code was created (early 90's), Borland did not supply a 32-bit Pascal which ran on the DOS and Win3.11 machines which most people used. To enjoy the 32-bit advantage, we had to use Borland-like Pascals, which were generally a few years behind Borland's development of the Pascal language. The SparSolv code was often run on Frontier Software 32-bit Pascal (TP5.0 compatible). There are/were a number of similar products (FPK, TMT, GNU, VP/OS2, StoneyBrook, etc).

For November 1999,a problem with a missing "Timer" unit has been fixed. A new,
Windows, program TESTMTX has been added which reads in and solves matrices in the
'MatrixMarket' MTX format. Three MTX files are supplied. TESTMTX allows you to compare
SparSolv with a similar unit, SparLin, obtainable from:

http://www-rab.larc.nasa.gov/nmp/fNMPhome.htm

Sparlin comes with several more MTX files not included with SparSolv.

Although considerably less sophisticated, SparSolv compares quite well with SparLin (which is a port of a NetLib C routine). The following table compares the times and fillins (non-zeroes added) of the two routines for several MTX files:

matrix size elements sparlin sparlin sparsolv diag option secs fills secs fills secs fills mahindas 1258 7687 0.6 (2292) 9.0(118815) 0.1 (3933) 1138_bus 1138 4059 0.3 (1489) 0.0 (1548) 0.0 (0) memplus 17758 126155 98.0 (25474) 4.2 (29259) 60.1 (19944) sherman2 1080 23099 31.2(207920) 46.9(326717) 3.8(137530) test1 2500 19992 222.5(700848) 263.8(950317) 45.0(746263) Notes: Times are on a PIII 350mhz (in 2007, using 3.2 GHz Pentium 4, times were 1/5 of above) SparLin has a "use diagonal" option which sometimes helps.

On the whole SparSolv seems faster: I like to think this is because it was developed in Pascal, while Sparlin was originally in C, a language that slows thinking. However, what we should really learn from comparisons such as this, is that no one sparse algorithm is best for all matrices. Sparlin seems to perform very well on diagonally dominant matrices (like Memplus), less well on more random matrices (like test1).

**September 2002:** minor revisions to accept compiler directives suitable for
Delphi 6. SparSolv seems to run fine in Delphis 3-6, either as part of a Windows program
or in a console program {$APPTYPE console}.

A third MTX example matrix is now included in the package.

Added example program ASOLVER.DPR contributed by Alex Jakushev, showing how to add
exception handling.

**August 2005:** minor revisions to accommodate Free Pascal Compiler version 2.0.0
[2005/05/08].

**January 2007:** compiler directives simplified -- they now assume that you are
using some 32-bit Delphi; if you are not, some hand editing will be needed.
[2007/19/01].

Complaints and suggestions to Mark.Horridge

email: Mark.Horridge@gmail.com

This page courtesy of my work-place: Centre of Policy Studies at Monash University

Last updated 6 Jul 2009.

The "Numerical Analysis in Pascal" website was maintained till 2005 by Mark Vaughan at at www-rab.larc.nasa.gov/nmp/fNMPhome.htm, but is no longer available there. You can download here a zip archive of the site. Unpack the zip into a new folder and browse fNMPhome.htm. Although many older external links are broken now, much useful information remains.