Re: Better quadrature routine in octave

[Jaroslav Hajek]

On Thu, Apr 1, 2010 at 1:36 PM, Pedro Gonnet <address@hidden> wrote:
> On Wed, 2010-03-31 at 14:09 +0200, Jaroslav Hajek wrote:
>> I advise you not to do it until it becomes apparent that the m-file
>> version can't be optimized to provide satisfactory performance.
>> Having algorithms in m-files has certain advantages, in particular
>> that it is easier for people to play with them.
>
> Ok, I was going to do it anyway for the GNU Scientific Library, but this
> somewhat takes the urgency out of it.
>
>> If you come up with a version that can replace quad (if I understood
>> correctly, it suffices to add handling of singular points)
>> and adapt it somewhat to the coding standards (at least the major
>> points, i.e. consistent comments, specific block ends), I can help you
>> with optimizations, should any be needed.
>
> That's exactly what I've done. I've also added the handling of +/-Inf in
> the interval, albeit in a hackish way, as well as some documentation and
> some test functions that I shamelessly stole from quadgk.
>
> I haven't done any extensive tests with the new code though! Will do
> this as soon as possible using the examples from the paper, so what I've
> attached here is NOT FINAL, but just a teaser to get some comments
> regarding form and style!
>
> Some preliminary observations: The code is much slower than it is in
> Matlab, probably due to the whole JIT thing. I can only guess as to
> where performance is lost -- I'm guessing a good part is lost in the
> "heap" which is not a heap.

I would not expect so. In Octave, structs gather objects using
(reference-counted) pointers, so manipulating struct arrays generally
does not involve moving much memory, even if the stored elements are
quite big.

> Back when I implemented it in Matlab, using
> and unsorted array of structs was faster than building and maintaining
> an explicit heap. Is there any way to profile the code in Octave to see
> if this is really where I should focus my attention?
>

No, there's no profiler yet (and I don't see it come for quite some time).

At a first glance I see a couple of vectorization possibilities,
though I can't be sure what the general impact will be.
Could you provide several test inputs (such as those you used) for
testing the speed?


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