Saturday, December 5, 2009

LS Solver and single precision

I was looking to add some facilities to the LS solver. In the process I was shown how the solver was built for inclusion in the GAMS distribution. Turns out parts of it are incorrectly compiled using single precision. I just about fell off my chair with astonishment when I saw this. After recovering somewhat, I investigated a bit what the damage could be. Luckily the main linear algebra routines are compiled correctly, so the main results (estimates) are fine. The affected pieces involve the p-values and related statistics. Comparing the results with the correctly compiled version from http://www.amsterdamoptimization.com/statistics.html, the differences seem minimal. Of course at this point it is better to use this version instead of the version included in the GAMS distribution. Let’s hope this will be fixed soon.

The background is some older code from TOMS (Transaction on Mathematical Software), algorithm TOMS708 for calculating the incomplete Beta Integral. This code was originally developed for the CDC 6600 architecture, which stores single precision numbers in 60 bits. For current IEEE machines that means this should be compiled as double precision. Similar issues always popped up for code developed on Cray super computers (the old vector machines had 64 bit single precision numbers). Luckily all fortran compilers have command line switches to make it easy to compile such code with double precision.

I encountered this brochure for the Control Data Cyber 170-750: the machine I used to work on (mainly Pascal and Fortran; the machine was running an operating system called NOS/BE, see http://en.wikipedia.org/wiki/CDC_Cyber).