[ublas] solving a linear system Ax=b with A singular

From: Marchesi Mauro (Marchesi_at_[hidden])
Date: 2007-06-01 09:10:12

Hello.

I tried to solve a linear system Ax=b by using the following code:

        ublas::matrix<double> A = ...;
        ublas::vector<double> b = ...;

        ublas::permutation_matrix<size_t> pm ( A.size1() );
        ublas::matrix<double> M = A;
        ublas::vector<double> x = b;
        ublas::lu_factorize( M, pm );
        ublas::lu_substitute( M, pm, x );

This works fine, unless A is singular.
However it turns out that, because of the way the matrix is built, the
rank of the coefficient matrix A equals the rank of the augmented matrix
Ab.

So, even in this case, the linear system admits solutions, but I can't
figure out how to get them by using the uBLAS library.

I'd like to get the set of solutions in any form: for example as a
vector x0 along with a set of vectors {x1,...,xN} such that every
solution of the system can be expressed in the form x0 + k1*x1 + ... +
kN*xN (clearly the set {x1,...,xN} is empty iff A is non-singular).

Can anyone help me?

Thanks in advance,
Mauro Marchesi