Re: [ublas] Error in ublas_gesv.cpp?

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