# 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 =

-1.0635e+013
-1.3903e-002
1.0635e+013
-2.1269e+013
3.1904e+013
-4.2539e+013
5.3173e+013
-6.3808e+013
7.4443e+013
-8.5077e+013
9.5712e+013
-1.0635e+014
1.1698e+014
-1.2762e+014
1.3825e+014
-1.4889e+014
1.5952e+014
-1.7011e+014
1.8012e+014
-1.8012e+014

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

-0.013903
0.996094
1.664063
3.410156
3.656250
5.179688
6.039063
6.906250
8.000000
9.234375
10.000000
11.375000
12.140625
12.562500
13.812500
14.531250
15.312500
15.250000
17.750000
19.000000

octave-3.0.3.exe:33>

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.

Regards,
Thomas