Seg-fault in gsl-0.5/integration/qpsrt.c

[Brett Viren <bviren at superk dot physics dot sunysb dot edu>]

Hello.

I guess this is for Brian Gough.

I think I spotted a problem in gsl-0.5/integration/qpsrt.c.

It turned up when using gsl_integration_qag to integrate a 5
dimensional integral containing various gaussian and poisson functions
along most dimensions.  

I had a bit of trouble getting gdb to halt just before failure because
the segfault comes after several thousand calls to
gsl_integration_qag.  Also, I had some problems understanding what was
going on due to inlined functions and optimization (can these be
turned off?).  But anyways, letting the segfault occur when running
under gdb pointed out the err as beeing at qpsrt.c:79 which is
right...

  while (errmin >= elist[order[k]] && k > i - 2)
    { 
      order[k+1] = order[k] ; 
      k-- ; 
    }			/* <---- ...here */

I don't understand why gdb sees the segfault at the end of the while
loop and I am also a bit confused at the logic in this file, but I
tried reversing the two tests, namely:

  while (k > i - 2 && errmin >= elist[order[k]]) 

and the segfault no longer occured.  To me, this seems the proper
order regardless, (keeping in mind I don't really understand the logic
of qpsrt.c <grin>).  There could also be the same problem at the while
loop just above on line 61:

  while (errmax < elist[order[i]] && i < top) 

But it was not being tickled by my code.

I tried to better isolate the problem to a test code, but was unable
to duplicate the condition.  If anyone wants my whole code to test it,
I'll supply it.


Thanks,
-Brett.

PS: I am not on this list, so I'll check for any responce in the
archive, or else just CC me.