Boost logo

Ublas :

Subject: Re: [ublas] Error in ublas_gesv.cpp?
From: Thomas Klimpel (Thomas.Klimpel_at_[hidden])
Date: 2011-03-27 12:19:30

Richard Stanton wrote:
> I just downloaded the latest boost numeric bindings from the boost sandbox,
> and tried to compile and run the sample program ublas_gesv.cpp, adding
> some lines to print out the solutions obtained to the original equation A x = b.

That's good, because I was expecting more or less that the tests might be used in that way.

> If I make the matrix small, e.g., size below 10, things seem to work fine
> – when I take the original A matrix and multiply by the proposed solution,
> I do indeed get the original b vector (note – I make copies of both A and b
> before testing this). However, at size 20, for example,
> here’s what I get for the double version of the code:
> x = [20](-1.36885e+14,-0.015625,1.36885e+14,-2.73771e+14,4.10656e+14,
> -5.47541e+14,6.84426e+14,-8.21312e+14,9.58197e+14,-1.09508e+15,
> 1.23197e+15,-1.36885e+15,1.50574e+15,-1.64262e+15,1.77951e+15,
> -1.91639e+15,2.05324e+15,-2.18962e+15,2.31842e+15,-2.31842e+15)
> B = [20](0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)
> Prod = [20](-0.015625,1,2,3.0625,4,5,6.125,7,7.625,9.5,10,10.25,
> 10.25,13,13.25,15.75,15,21.5,11.5,19.5)
> The last line should just repeat the line above, if everything went according to plan.

Thanks for reporting this. I should probably fix this test sooner or later (use a well conditioned matrix, add checks for correctness of result). What happens here is that the matrix has a very bad condition number. See below how I tried to with octave to handle this matrix:

octave-3.0.3.exe:27> a = diag([0:19]);
octave-3.0.3.exe:28> a(2:20,1:19) += diag([1:19]);
octave-3.0.3.exe:29> a(1:19,2:20) += diag(ones(19,1));
octave-3.0.3.exe:30> b = [0:19]';
octave-3.0.3.exe:31> x = a\b
warning: matrix singular to machine precision, rcond = 3.98706e-019
warning: attempting to find minimum norm solution
x =


octave-3.0.3.exe:32> a*x
ans =



As a final note, don't be disappointed that the test won't be immediately fixed. I will take care of it when I have time.