[SciPy-Dev] Quadpack conversion to Fortran

Pauli Virtanen pav@iki...
Sun Sep 8 14:53:24 CDT 2013

Hash: SHA1

08.09.2013 07:43, Brian Lee Newsom kirjoitti:
> My name is Brian Newsom and I am just beginning to use/work on the
> scipy library, specifically the quadpack within the integration
> pack. Because of the recursion and wrapping methods used now,
> running nquad (scipy 0.13) function is slow-the fortran and python
> are nested in a way that requires switching back and forth between
> the two often.  Ideally this inefficiency could be overcome by
> translating the main interface (quadpack.py) into fortran and just
> having a minimal python wrapper that is not used in the recursion.

The first question: did you measure the overhead involved in the
Python code, and determined that it dominates the computation time? If
not, that would be the first thing to do, otherwise it's possible to
end up with a wild goose chase.

On the face of it, it does not sound likely to me that much will be
gained by writing the recursion in a lower-level language. Suppose
that each level requires N >> 1 evaluations of the function. The
innermost function will be called in total N^3 times, the next level
N^2 times and the outer level N times. The inner level involves N^3
calls back to Python, which are not removed by writing the recursion
in Fortran. Since N^3 >> N^2, the innermost callback overhead
dominates the run time. Of course, the N^2 term has a bigger prefactor
than N^3, so this argument is not fully convincing, but it would be
best to measure first.

IIRC, the `quad` integration routine can directly call ctypes function
pointers, skipping Python overhead. (See the Python ctypes module
documentation for details on how to define them.) This requires Scipy
> = 0.11.0, however.

There have been plans to improve the interoperability further (with
Cython, f2py, and other ways) to interact directly with low-level code
in Scipy's computational routines. Putting effort into this direction
and trying to speed up the innermost part would seem to me the most
sensible direction of further work.

Of course, the above is written assuming your benchmarks don't show me

I don't believe anyone is currently working on `nquad`, so go ahead.
However, note that due to the current situation with Mingw Gfortran
compiler on Windows, there may be issues in accepting Fortran 90 code.

Best regards,

- -- 
Pauli Virtanen
Version: GnuPG v1.4.12 (GNU/Linux)


More information about the SciPy-Dev mailing list