So, contrary to all good opinion and intention, I find myself back in the guts of the 10 year old C / C++ monstrosity that is BLAST. For those who aren’t bioinformatics geeks, BLAST is the Basic Local Alignment Search Tool. It’s the very most popular piece of software ever created for analyzing genetic sequence data. It’s also a big pile of C, written by biologists for biologists. One of those pieces of software inside of which you really don’t want to look.
So anyway, we have one architecture (Alpha E40, running Redhat 7.2, thanks for asking) for which BLAST barfs on a particular alignment. If I feed it two very specific sequences it dumps core. So far, I’ve rebuilt without optimization, upgraded gcc to the very most recent version, and examined the core file. This showed me that it’s an “Arithmetic error” in a particular function. Turns out it’s the line:
return searchsp * exp((Nlm_FloatHi)(-Lambda * S) + kbp->logK);
Everything is cool until that outermost multiplication by a very small value. Example idiot prints:
Entering BlastKarlinStoE_simple
^^ searchsp: 351649.000000 Lambda: 1.374063 S: 489 logK: -0.341642
^^ -Lambda * S: -671.916867
^^ exp((-Lambda * S) + kbp->logK): 1.101132e-292
^^ searchsp * the above: 3.872119e-287
Leaving BlastKarlinStoE_simple
Entering BlastKarlinStoE_simple
^^ searchsp: 351649.000000 Lambda: 1.374063 S: 518 logK: -0.341642
^^ -Lambda * S: -711.764697
^^ exp((-Lambda * S) + kbp->logK): 5.446862e-310
Floating exception (core dumped)
The part that gets me is that this works *fine* on other architectures. Unless I get this done in pretty short order, it would have been cheaper for me to upgrade the system on which the code runs than to spend my time debugging it.
Hate.
