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:
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:
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.
Unit SparSolv also makes available 3 LongInt variables to assist in monitoring and controlling memory use:
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!
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:
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].
September 2020: I Googled "sparsolv", and learnt that 2 or 3 people had adapted SparSolv to run with later Delphi versions. Added note above on complex-number linear systems.
Complaints and suggestions to Mark.Horridge
email: Mark.Horridge@gmail.com
This page courtesy of my work-place: Centre of Policy Studies
Last updated September 2020.
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.