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> // 1
#include <rw/math/genfact.h>
int main()
{
RWGenMat<float> A; // 2
RWMathVec<float> b;
std::cin >> A >> b; // 3
RWMathVec<float> x = solve(A,b); // 4
std::cout << x << std::endl;
return 0;
}
Here is a line-by-line description of what this program does:
2x2
[ 4 5
1 3 ]
[ 2 1 ]
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!
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 (LU.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.