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,
> heres what I get for the double version of the code:
> x = (-1.36885e+14,-0.015625,1.36885e+14,-2.73771e+14,4.10656e+14,
> B = (0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)
> Prod = (-0.015625,1,2,3.0625,4,5,6.125,7,7.625,9.5,10,10.25,
> 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
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.