float vs double .... what a difference

[Klaas]

Hi,

Using lapack 3.1 fortran version linked into a c++ program I tried to test something simple and was amazed by the result.

when inverting a symetrical 4x4 real matrix by means of spptrf and spptri or dpptrf and dpptri the result is so much different.

These are the elements of the matrix i tried to invert
a(0,0) = 1.25331
a(0,1) = 1.15695
a(1,1) = 1.25331
a(0,2) = 1.15695
a(1,2) = 0.910093
a(2,2) = 1.25331
a(0,3) = 1.19817
a(1,3) = 0.980975
a(2,3) = 1.24706
a(3,3) = 1.25331

when i stored a in single precision and inverted the the matrix by means of spptrf and spptri the result was

Code: Select all

 2338.06   -559.483  3647.85  -5426.94
 -559.483   141.22   -834.062  1254.24
 3647.85   -834.062  5977.17  -8781.9
-5426.94   1254.24  -8781.9   12945.4


when i did everything in double precision i obtained

Code: Select all

 2339.55   -559.834  3650.2   -5430.43
 -559.834   141.303  -834.619  1255.06
 3650.2    -834.619  5980.89  -8787.41
-5430.43   1255.06  -8787.41  12953.6

To be honest ... WOW! what a bunch of bits can do.

[sven]

Your matrix has a condition number with respect to inversion of approximately 10^5, so you can expect to lose about 5 digits of accuracy in the computed solution. Your results seem to show that nicely.

Best wishes,

Sven Hammarling.

[Julien Langou]

Sven's right, your dealing with an ill-conditioned problem. So your result is not surprising.

However, I want to make a remark.

Which solution is the best: the one from the single precision solve or the one from the
double precision solve? Be careful with the fast conclusion to say that double precision is
better than single precision. You can draw this conclusion only if you can ensure to know
your input data (i.e. matrix A) more precisely than the single precision accuracy.

For example, assume that a(0,0) is not

Code: Select all

a(0,0) = 1.2533100,

but a(0,0) is

Code: Select all

a(0,0) = 1.2533101.

Solve both in double precision arithmetic, and well, you'll have the same kind of variation
than between single and double. In other words, knowing your input data to 8 digit
accuracy is only enough to determine the three first digits of the inverse of your matrix
(since the condition number is 10^5). In other words, the single precision answer has the
same quality as the double precision answer if your input data is only known to 8 digit
accuracy.

-j