Solving Several Systems of Equations

SourcePro : Linear Algebra Module User’s Guide : Factorizations : 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> // 1

#include <rw/lapack/pdfct.h>

int main()

{

RWSymMat<double> A; // 2

std::cin >> A;

RWPDFact<double> fact(A); // 3

if (fact.good()) { // 4

RWMathVec<double> b; // 5

while(std::cin >> b) { // 6

RWMathVec<double> x = solve(fact, b); // 7

std::cout << x << std::endl;

}

return 0;

}

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.