A Simple Simultaneous Equation Solver
by K.B. Williams
Introduction
This article describes a function that performs the following tasks:
1. Inverts a matrix A (say).
2. Solves the accompanying system of right-hand sides, if
any.
3. Computes the value of the determinant of A.
What differentiates this program from the run-of-the-mill equation solvers are the following features:
1. The indexes of matrix elements are
computed efficiently — a two-dimensional array is treated
as singly dimensioned.
2. Any number of right-hand sides can be included.
3. The matrix may be embedded in a larger array.
This third feature has proved invaluable to me over the years. I have used it to solve pieces of a system of equations. I have also used it in an algorithm that inverts a large matrix by the method of partitioning [1].
The name of the function is augvrt (Listing 1) . It started its journey down the path of scientific computation in the mid 1960s. The source was originally written in Fortran II and happened my way while I was laboring mightily in computer programming in Melbourne, FL, in 1966. I have carted this gem from Florida to Hawaii, then to California and finally to Oklahoma. Along the way it has worked its magic in the following languages:
1. Fortran II
2. Fortran IV
3. CDC Assembler (COMPASS)
4. Ratfor
5. ANSI Fortran
6. C
The method is straightforward good old Gaussian elimination [1] . You can find a description of this method in nearly every numerical analysis book ever written. The value of the determinant is calculated in the traditional way except that if underflow occurs a value of 1.0e-30 is stored.
One drawback that you might want to be aware of is that no row or column pivoting, partial or otherwise, is performed. I found no need for features such as these even when working on problems in reentry missile-data reduction and structural analysis.
Using augvrt
A useful utility must work as advertised and must be dead easy to use. Function augvrt is such a utility. The calling sequence is simple:
double augvrt(double *a, int nr, int nc, int nd);
where:
a is the address of matrix to be inverted
nr is the number of rows in matrix a
nc is the number of columns in matrix a (nc >= nr)
nd is the number of columns in the matrix that holds matrix
a (nd >= nc)
(The matrix to be inverted can be embedded in another.) If x = A - inv * b is being solved, the nc - nr augmenting vectors will be replaced by nc - nr solution vectors.
The quantity returned is the value of the determinant of matrix a. If zero is returned, the matrix was not inverted and no solutions were produced. If underflow occurs at any step in the computation of the value of the determinant, a default value of 1.0e-30 will be stored.
Test Cases
I have provided three test cases for your pleasure (on this month's code disk).
tstvrt.c — This test inverts a 4 x 4 matrix and produces a solution vector by multiplying the inverse by an augmenting vector held separately. The inverse and the solution vector are printed followed by the known results. You are expected to do your part by comparing the two.
tstvrt2.c — This test operates on a 4 x 5 matrix. The fifth column is the augmenting vector. The matrix is inverted and the solution vector in column 5 is produced simultaneously. The 4 x 5 result matrix is printed, followed by the known results.
tstvrt3.c — This test operates on a 4 x 6 matrix. The fifth column is the augmenting vector and the sixth is to be ignored. The matrix is inverted, the solution vector in column 5 is produced simultaneously, and column 6 is ignored. Output from augvrt is printed followed by the known results.
All test cases will let you enter a fudge factor to be added to the first element of the input matrix. This lets you experiment with augvrt to see how numerical problems can occur.
You can compile this suite of tests, once you have the source in place, by entering any of nmake, bmake, make, or whatever the name of your favorite make utility is. The makefile is set up to use Microsoft C/C++, but a simple change will let you do your thing with Borland C, or other fine compilers.
A separate manual page accompanies the code on this month's code disk. The test program requires several other functions, one of which is matmpy [2] . All functions are included on the code disk.
Acknowledgments
Lost in antiquity are the origins of the original code for augvrt and the test matrix used in the test programs. It has been 29 years since I first laid eyes on augvrt while working at D. Brown Associates in Melbourne, Florida. The test matrix has been with me for well over 20 years; it is buried in a numerical analysis book somewhere. Thanks to those unknown sources.
References
1. D.K. Faddeev and V.N. Faddeeva. Computational Methods of Linear Algebra (W.H. Freeman & Co., 1963).
2. K.B. Williams, "Magic Matrix Multiply," The C/C++ Users Journal, June 1995, pp. 19-24.
K.B. Williams received a B.A. in Mathematics from
California State University of Sacramento in 1955. He is a
veteran of forty years in the scientific applications
programming trenches and is the founder of King Baker
Enterprises, of Stillwater, Oklahoma. He consults in applied
numeric techniques. He can be reach at 802 South Ridge
Drive, Stillwater, OK 74074, or via e-mail at [email protected].
