Re: Better quadrature routine in octave

[Pedro Gonnet]

Hi again,

I just finished re-coding the quadrature routine in C for inclusion in
the GNU Scientific Library. The code is blazingly fast and passes all
the tests in the paper. It produces essentially the same results as the
Octave code, differing only in how it selects the interval with the
smallest error for removal: while Octave uses min(..), the C version
truncates the heap, giving slightly different (yet correct) results for
the more complex integrals in the test series.

If you've got GSL installed anywhere, you can compile the code with

         gcc -g -lm -lgsl -lblas cquad_test.c

The main routine in cquad_test.c implements the battery tests from the
paper.

The code uses the GSL only for its representation of the integrand and
for the error codes. These parts would have to be re-written anyway to
make it conform to the Octave C-language interface.

I've only just started reading the documentation on how to do this, e.g.
add C-language code to Octave. Any help would be greatly appreciated!

Cheers,
Pedro


On Thu, 2010-04-15 at 13:30 +0200, Jaroslav Hajek wrote:
> On Tue, Apr 13, 2010 at 4:54 PM, Pedro Gonnet <address@hidden> wrote:
> >
> > Hello again,
> >
> > Ok, I've checked Shampine's original paper on quadgk/quadva and there
> > are no real, i.e. systematic or parameterized, tests for the interval
> > transformation, so I wouldn't know what else to test my code against.
> >
> > How do you guys want to proceed on this? If anything needs
> > fixing/improving, just tell me what and I'll get to it.
> >
> > Cheers, Pedro
> >
> 
> I finally got into testing this, and for some functions of the test
> suite (e.g. floor(exp(x)) on 0..3) it is indeed a bit slow.
> I'm sure it can be improved, but I would start with simplifying the code.
> 
> Your code has a lot of repeated patterns. The quadruples of analogical
> vars like V_1, V_2, ... can be replaced by cell arrays, repeated
> similar code blocks can be moved into subfunctions. I started with
> some refactoring, the result is attached (this is *not* a working
> version) if you wish to continue.
> 
> Regarding the optimizations, there are some immediate targets. Stuff like
> 
>         % check for nans and clean them up
>         nans = [];
>         for i=1:length(fx), if ( ~isfinite (fx(i)) ), nans = [ i ,
> nans ]; fx(i) = 0.0; endif; endfor;
> 
>         ...
> 
>         % re-instate the nans
>         for i=nans, fx(i) = NaN; endfor;
> 
> is *much* better written in a vectorized way as
> 
>         % check for nans and clean them up
>         nans = find (! finite (fx));
>         fx(nans) = 0;
> 
>         % re-instate the nans
>         fx(nans) = NaN;
> 
> the repeated nested references are another potential source of
> sluggishness, I think lots of them can be eliminated.
> In general it's better to first simplify the code if possible, then
> start with optimizations.
> 
> regards
>

cquad.c
Description: Text Data

cquad.h
Description: Text Data

cquad_consts.h
Description: Text Data

cquad_test.c
Description: Text Data