[Numpy-discussion] Cython-based OpenMP-accelerated quartic polynomial solver

Chris Barker - NOAA Federal chris.barker at noaa.gov
Tue Sep 29 20:48:58 EDT 2015


This sounds pretty cool -- and I've had a use case. So it would be
nice to get into Numpy.

But: I doubt we'd want OpenMP dependence in Numpy itself.

But maybe a pure Cython non-MP version?

Are we putting Cuthon in Numpy now? I lost track.

-CHB

Sent from my iPhone

> On Sep 29, 2015, at 7:35 AM, Juha Jeronen <juha.jeronen at jyu.fi> wrote:
>
> Hi all,
>
> I recently developed a Cython-based, OpenMP-accelerated quartic (and cubic, quadratic) polynomial solver to address a personal research need for quickly solving a very large number of independent low-degree polynomial equations for both real and complex coefficients.
>
> For example, on an Intel i7-4710MQ, running with 8 threads this code solves approximately 1.2e7 quartic polynomial equations per second. (With OMP_NUM_THREADS=4 the solver gives approximately 8.9e6 solutions per second, so it seems HyperThreading in this case helps about 30%.)
>
> The algorithms for cubics and quadratics come from Numerical Recipes (3rd ed.), and the quartic problem is internally reduced to a cubic and two quadratics, using well-known standard tricks.
>
>
> Since to my understanding, for solving polynomial equations NumPy only provides roots(), which works one problem at a time, I'm writing to inquire whether there would be interest to include this functionality to NumPy, if I contribute the code (and clean it up for integration)?
>
>
> I have some previous open source contribution experience. I have authored the IVTC and Phosphor deinterlacers for VLC, and a modular postprocessing filter framework for the Panda3D game engine (posted at the forums on panda3d.org, projected for version 1.10).
>
> Currently the polynomial solver is available in a git repository hosted by our university, the relevant files being:
>
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve2.pyx
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve.py
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve2example.py
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/setup.py
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/compile.sh
>
> I'm working on a research grant, which allows me to hold the copyright; license is BSD.
>
> Polysolve2 is the fast Cython+OpenMP version, while polysolve is an earlier (slower) experimental version using only NumPy array operations. The example module contains a usage example, and setup (in the standard manner) instructs distutils and Cython as to how to build polysolve2 (including setting the relevant math flags and enabling OpenMP).
>
> (The setup script includes also some stuff specific to my solver for the Ziegler problem; basically, the "tworods" module can be ignored. I apologize for the inconvenience.)
>
>
> Thoughts?
>
>
> Best regards,
>
> -J
>
> -------------------------------------------------
> Juha Jeronen, Ph.D.
> juha.jeronen at jyu.fi
> University of Jyväskylä
> Department of Mathematical Information Technology
> -------------------------------------------------
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion



More information about the NumPy-Discussion mailing list