Linear Equation Solver in C++

• Download source code - 31.7 KB

Introduction

The program SolveLinearEquations solves simultaneous equations. The program takes a text file that contains equations for input and outputs the solution. This program was written years before I wrote the Linear Equations Solver in C# here. Other than that, this program has no graphical user interface and a slightly modified equation format, this program is very similar to the C# program.

The program uses the SparseArray template class to implement vectors and a matrix. The matrix uses the DoubleIndex class, which takes two integer indices, to implement a single key to use with the SparseArray class. The stored type for the matrix SparseArray is a double precision number. Each vector is just an instance of the SparseArray class that stores a double precision value with an integer index for a key.

The code uses a crude parser to parse the equations in the input file. It's crude because it requires a somewhat rigid input format and doesn't support parenthesis or math functions.

The input file is a text file that contains equations. There can be no more than 1024 characters on a line. An equation is terminated with a semicolon. If there is no semicolon at the end of a line, then the equation is continued on the next line.

The equations have the following format, which allows using the addition and subtraction sign to combine terms of the form:

[number][variable]

There must be a single equals sign in each equation.

Either the number or variable are optional and as many terms as desired can be combined using either the plus sign or a minus sign. The number can contain a decimal point and an exponent. Variables can
only contain alphabetic characters or the underscore character.

Floating point exponents are preceded by the ^ character instead of the usual E character to avoid any ambiguity with variable names. The following line shows a floating point number equal to 2.3 million.

X = 2.3^6

An equation must contain an equals sign.

An example set of equations is:

3  X + 4 Y = -5 Z;
X + Z = 10 Y;
X + Z = 42.5;

Note the blank line that separates the equations.

The program produces the solution for those equations:

X = 114.75
Y = 4.25
Z = -72.25

Another example set of three equations might be:

MARYS_AGE=2BOBS_AGE ; BOBS_AGE = 3 CATHYS_AGE;

CATHYS_AGE = 4;

or another example:

HEIGHT = 5 + 10; 

If You Forget The Equation Format

Running the program with the following command line will give the same information as above regarding equation format required in the input file.

SolveLinearEquations -h

Background

When I was in college, I wrote a circuit analysis program in Fortran. I needed a way to solve simultaneous equations, and I stumbled on the following book and algorithm:

"The solution of ill-conditioned linear equations", by J. H. Wilkinson, "Mathematical Methods For Digital Computers" Volume 2, Edited by Anthony Ralston and Herbert S. Wilf, 1967, John Wiley and Sons, pages 65-93.

An ill-conditioned set of equations is a set of equations that is either difficult or impossible to solve using the given floating point precision. Equations are ill-conditioned when two or more equations define almost parallel lines, or in the case of more than two dimensions, almost parallel planes or almost parallel hyperplanes. An example of ill-conditioned equations would be:

3 10^-12 X + Y = 0.7

X + Y = 0.9

I was lucky to stumble on this particular reference. While this book gives a standard implementation of Gaussian Elimination with partial pivoting, as is taught in Computer Science Linear Algebra courses, this algorithm also determines whether an accurate solution can be found. This is done using both the matrix norms and a constant that is set based on the number of bits in the mantissa of a floating point number.

I've rewritten this algorithm several times. I wrote this in Fortran, C++ using simple arrays, again in C++ using a sparse container classes, and finally in C#.

Back in the late 1970s, I ran this algorithm on a DEC-10 computer and solved 1000 equations with 1000 variables in about 30 seconds. Today, with that original C code, a problem that size on a PC runs in the blink of an eye. This code is slower because it uses a SparseArray, but it is still very fast.

File List

• SolveLinearEquations.cpp - The main program
• LinearEquationParser.cpp - Parses the input file that contains equations
• LinearEquationParser.h - The parser header file
• CharString.cpp - A string class
• CharString.h - The string class header file
• MatrixPackage.cpp - Contains the linear equations solver function
• MatrixPackage.h - The linear equation solver header file
• SparseArray.h - The SparseArray template file
• DoubleIndex.cpp - A class to make two indices into one key
• DoubleIndex.h - The header file for DoubleIndex
• SolveLinearEquations.vcproj - A Visual Studio 2008 project file
• SolveLinearEquations.sln - A Visual Studio 2008 solution file

About the SolveLinearEquations Function

A set of Linear equations are represented by the matrix equation:

aMatrix xVector = bVector 

The aMatrix and bVector are given, and the xVector is the solution.

The first example set of equations given above can be rewritten as:

3 X +  4 Y + 5 Z = 0
1 X - 10 Y + 1 Z = 0
1 X +  0 Y + 1 Z = 42.5   

The matrix form of these equations is:

| 3   4  5 | | X |   |   0  |
| 1  10  1 | | Y | = |   0  |
| 1   0  1 | | Z |   | 42.5 |

The aMatrix is the square matrix on the left. The bVector is on the right.

The xVector, which contains the variable names, which are the unknowns, is in the middle.

To solve these equations, call the SolveLinearEquations function in the MatrixPackage namespace. This function is defined in file MatrixPackage.h and implemented in file MatrixPackage.cpp.

Status_T SolveLinearEquations(unsigned int number_of_equations,
                              const SparseMatrix & a_matrix,
                              const SparseVector & b_vector,
                              SparseVector & x_vector);

The intent of naming those files MatricPackage was with the intent of one day implementing a larger set of matrix operations. All that I needed at the time was an equation solver, and that, and some error code, is all that is in the matrix package now.

The xVector, which will store the solution, is the last argument. The number of equations is the size of one dimension of the square matrix.

The program will also indicate if a set of equations is 'singular' to working accuracy. A singular set of equations has no single solution because two or more equations are merely a multiple of the other equation, such as:

X + Y = 7
2X + 2Y = 36

Even if the second equation was, "2X + 2Y = 14", so it was consistent with the first equation, there is no single solution to the two equations, and the program would report that the equations were singular to working accuracy.

An Important Parameter For Determining If Equations Are Ill-Conditioned

The file MatrixPackage.cpp contains a constant that is set based on the number of bits in the mantissa of a double precision floating point number.

const double f_SMALL_FLOAT = 5.69E-14;

If you port the linear equation solver to another platform, then be sure to adjust this constant. The code contains the following comment regarding this constant.

// The original implementation was on a computer where the floating
// point mantissa was 39 bits.  For that system, the value below
// was set to 2.92E-11, which is just slightly larger than
// 1/(2^35), which is 2.91E-11.  For my Intel system, the mantissa
// of a double precision floating point number is 48 bits, so
// the value is set to slightly greater than 1/(2^44).  1/(s^44)
// evaluates to 5.68E-14, so the value 5.69E-14 is used here.

Why Sparse Containers are used for the Vectors and Matrices.

In the late 1970s, I implemented this algorithm in Fortran for a circuit analysis program. The arrays were hard-coded sizes. So, to solve 1000 by 1000 system of equations requires two matrices that both had 1000² entries, or a million entries! Years later, I recoded this in the C language, again with hard-coded sizes.

Most real-world problems either require matrix sizes much smaller than 1000 by 1000, or the matrix is sparse. The circuit analysis program typically required roughly 5 terms in each equation, where each term corresponds to a column of the "a" matrix. So, a 1000 by 1000 matrix for such a problem would only have 1000 by 5, or 5000 non-zero entries. So, a sparse container would store only 5000 double precision values, whereas a full 1000 by 1000 matrix would require a million double precision values, or 200 times as many values!

For very large problems, the enormous space savings is worth the degradation in run-time performance.

Also, I realized that it's pretty easy to take an algorithm implemented with sparse containers and convert it to used fixed-size arrays, whereas doing the reverse is a lot of work. So, someone could convert this to use a lot of space, and possibly be faster. With modern processor caching, depending on bus-speeds and processor architecture, making an algorithm use a lot of memory can make the algorithm slower, even when the number of instructions executed is made smaller.

Using the Code

Create a file that contains only equations in the format mentioned earlier in this article, or use the provided file equations.txt. Enter:

SolveLinearEquations equations.txt 

The program will output an alphabetically sorted list of variable names and values. For file equations.txt, the output will be:

Linear Equation Solver - Version 2.01
Copyright (C) William Hallahan 2001-2013.
Ann = 2
Joe = 8
Mary = 12
Rita = 18
Tom = 70

Points of Interest

Only ASCII builds are supported on Windows. The CharString class does support UNICODE, but most of the other code would have to be ported. The most difficult code to port would be the code that uses ofstream for the input file. This is possible though, but it requires more work to make it portable across platforms.

I have not tested the code in SparseArray.h and MatrixPackage.cpp (.h) on Linux using the g++ compiler. I expect that if the code doesn't already compile and work on Linux, the changes required to make it work there would be minimal and simple.

History

• Initial post