4.3 Solving Several Systems of Equations

In this section, a slightly more complicated example program is given. This program reads a matrix A from the input stream, continues to read in right-side vectors, b, and solves the system of equations Ax=b until an end of file is reached. To make it more complex, we assume that the matrix A is positive definite symmetric. Such a matrix can be constructed for test purposes, if you'd like to try this program, by choosing a nonsingular matrix X and forming A=XTX.

#include <iostream.h>                             // 1
#include <rw/dpdfct.h>

main()
{
    DoubleSymMat A;                               // 2
    cin >> A;
    DoublePDFact fact(A);                         // 3
    if (fact.good()) {                            // 4
     RWMathVec<double> b;                         // 5
     while(cin >> b) {                            // 6
     RWMathVec<double> x = solve(fact, b);        // 7
     cout << x << endl;
     }
    }
}

Here is a description of the main points in the program:

//1	The header files for the input/output stream classes and the positive-definite factorization class are included.
//2	The double precision symmetric matrix `A` is defined and then read in from standard input.
//3	Here a positive-definite factorization of the matrix `A` is constructed and given the name `fact`. For this to work, the matrix `A` must be, in fact, positive-definite. If not, the results (`fact`) are invalid.
//4	Here we test to see if all went well. If so, then...
//5	We define a vector `b`...
//6	We keep reading vectors from standard input until an end-of-file is encountered, and then...
//7	We use the factorization to solve each set of equations `Ax=b`.

This program can be written without constructing the factorization explicitly by replacing //7 with:

RWMathVec<double> x = solve(A, b);

This is not advisable because the factorization must be computed every time through the loop even though it is not changing. Because constructing the factorization is typically the most expensive part of solving a set of equations, the result is a much slower program.