Re: Better quadrature routine in octave

[Jaroslav Hajek]

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

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

cquad.m
Description: Text Data