4.2 Solving a System of Equations

The most basic problem of linear algebra is to solve a system of equations. Here is a simple program to do this, using float precision and general square matrices:

#include <iostream.h>                     // 1
#include <rw/math/genfact.h>

main()
{
    RWGenMat<float> A;                    // 2
    RWMathVec<float> b;
    cin >> A >> b;                        // 3
    RWMathVec<float> x = solve(A,b);      // 4
    cout << x << endl; 
}

Here is a line-by-line description of what this program does:

//1	The first `#include` statement allows the use of iostreams. The second `#include` statement includes the header file `genfact.h`, which declares the float general factorization type *RWGenFact<float>*.
//2	Here the variable `A` is defined as a general floating point matrix, `RWGenMat<float>`, and `b` is defined as a floating point vector, `RWMathVec<float>`.
//3	The variables `A` and `b` are read in from the standard input. Here is a possible input file showing the format: 2x2 [ 4 5 1 3 ] [ 2 1 ] In this case, the 2 x 2 matrix: and the two-element long vector (2, 1) are read in.
//4	This is the key statement in the program. The vector `x` which satisfies the equation `Ax=b` is returned by the `solve()` function. We will discuss how this works in more detail later.

At this point, you notice that the factorization class, RWGenFact<float>, is not used in the program! Or is it? If you look in the header file genfact.h, you can see that the solve function is declared to take an argument of type const RWGenFact<float>&. The C++ compiler constructed a temporary variable of type RWGenFact<float> for you. This temporary variable was constructed using the RWGenFact<float>(const RWGenMat<float>&) constructor and then passed to the solve function. Class RWGenFact<float> was used after all!

4.2.1 Error Checking

At this point, you may ask: "What happens in the example program if the matrix A is singular?" The answer is less than satisfactory: a runtime error occurs. Fortunately, such errors are easily prevented by testing the factorization. This is done by replacing //4 with the code:

RWGenFact<float> LU(A);
if (A.fail()) { /* Your error handler goes here*/ }
RWMathVec<float> x = solve(LU,b);

The fail() member function exists for all the factorization types. If it returns TRUE, calling solve() will likely result in a runtime error. In this example, you could also use the member function isSingular() to test for an error condition, but some factorizations may fail even if the matrix is not singular, so calling fail() is a better method of checking.

In addition to singularity, other situations can induce runtime error; for example, attempting to construct a positive-definite factorization from a matrix which is not positive-definite.