> On Sat, Mar 14, 2009 at 3:58 PM, Sturla Molden <stu...@molden.no> wrote:
> There is also a ticket (#579) to add an implementation of the Bluestein
> algorithm for doing prime order fft's. This could also be used for zoom
> type fft's. There is lots of fft stuff to be done. I wonder if some of it
> shouldn't go in Scipy? I think David added some dcts to Scipy.
I am not changing or adding algorithms for now. This is just to prevent
NumPy from locking up the interpreter while doing FFTs.
The loops that are worth multithreading are done in C in
fftpack_litemodule.c, not in Python in fftpack.py. I have added OpenMP
pragmas around them. When NumPy gets a build process that supports OpenMP,
they will execute in parallel. On GCC 4.4 is means compiling with -fopenmp
and linking -lgomp -lpthread (that goes for mingw/cygwin as well).
The init function seems to be thread safe. cffti and rffti work on arrays
created in the callers (fftpack_cffti and fftpack_rffti), no global
objects are touched.
I'm attaching a version of fftpack_litemodule.c that fixes most of what I
mentioned.
Sturla Molden
#include "fftpack.h"
#include "Python.h"
#include "numpy/arrayobject.h"
static PyObject *ErrorObject;
/* ----------------------------------------------------- */
static char fftpack_cfftf__doc__[] = "";
PyObject *
fftpack_cfftf(PyObject *NPY_UNUSED(self), PyObject *args)
{
PyObject *op1, *op2;
PyArrayObject *data;
PyArray_Descr *descr;
double *wsave, *dptr;
npy_intp nsave;
int npts, nrepeats, i;
if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
return NULL;
}
data = (PyArrayObject *)PyArray_CopyFromObject(op1,
PyArray_CDOUBLE, 1, 0);
if (data == NULL) {
return NULL;
}
descr = PyArray_DescrFromType(PyArray_DOUBLE);
if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
goto fail;
}
if (data == NULL) {
goto fail;
}
npts = data->dimensions[data->nd - 1];
if (nsave != npts*4 + 15) {
PyErr_SetString(ErrorObject, "invalid work array for fft size");
goto fail;
}
nrepeats = PyArray_SIZE(data)/npts;
NPY_SIGINT_ON;
Py_BEGIN_ALLOW_THREADS
#pragma omp parallel for private(i,dptr) default(shared)
for (i = 0; i < nrepeats; i++) {
dptr = (double *)data->data + i*npts*2;
cfftf(npts, dptr, wsave);
}
Py_END_ALLOW_THREADS
NPY_SIGINT_OFF;
PyArray_Free(op2, (char *)wsave);
return (PyObject *)data;
fail:
PyArray_Free(op2, (char *)wsave);
Py_DECREF(data);
return NULL;
}
static char fftpack_cfftb__doc__[] = "";
PyObject *
fftpack_cfftb(PyObject *NPY_UNUSED(self), PyObject *args)
{
PyObject *op1, *op2;
PyArrayObject *data;
PyArray_Descr *descr;
double *wsave, *dptr;
npy_intp nsave;
int npts, nrepeats, i;
if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
return NULL;
}
data = (PyArrayObject *)PyArray_CopyFromObject(op1,
PyArray_CDOUBLE, 1, 0);
if (data == NULL) {
return NULL;
}
descr = PyArray_DescrFromType(PyArray_DOUBLE);
if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
goto fail;
}
if (data == NULL) {
goto fail;
}
npts = data->dimensions[data->nd - 1];
if (nsave != npts*4 + 15) {
PyErr_SetString(ErrorObject, "invalid work array for fft size");
goto fail;
}
nrepeats = PyArray_SIZE(data)/npts;
dptr = (double *)data->data;
NPY_SIGINT_ON;
Py_BEGIN_ALLOW_THREADS
#pragma omp parallel for private(i,dptr) default(shared)
for (i = 0; i < nrepeats; i++) {
dptr = (double *)data->data + i*npts*2;
cfftb(npts, dptr, wsave);
}
Py_END_ALLOW_THREADS
NPY_SIGINT_OFF;
PyArray_Free(op2, (char *)wsave);
return (PyObject *)data;
fail:
PyArray_Free(op2, (char *)wsave);
Py_DECREF(data);
return NULL;
}
static char fftpack_cffti__doc__[] ="";
static PyObject *
fftpack_cffti(PyObject *NPY_UNUSED(self), PyObject *args)
{
PyArrayObject *op;
npy_intp dim;
long n;
if (!PyArg_ParseTuple(args, "l", &n)) {
return NULL;
}
/*Magic size needed by cffti*/
dim = 4*n + 15;
/*Create a 1 dimensional array of dimensions of type double*/
op = (PyArrayObject *)PyArray_SimpleNew(1, &dim, PyArray_DOUBLE);
if (op == NULL) {
return NULL;
}
NPY_SIGINT_ON;
Py_BEGIN_ALLOW_THREADS
cffti(n, (double *)((PyArrayObject*)op)->data);
Py_END_ALLOW_THREADS
NPY_SIGINT_OFF;
return (PyObject *)op;
}
static char fftpack_rfftf__doc__[] ="";
PyObject *
fftpack_rfftf(PyObject *NPY_UNUSED(self), PyObject *args)
{
PyObject *op1, *op2;
PyArrayObject *data, *ret;
PyArray_Descr *descr;
double *wsave, *dptr, *rptr;
npy_intp nsave;
int npts, nrepeats, i, rstep;
if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
return NULL;
}
data = (PyArrayObject *)PyArray_ContiguousFromObject(op1,
PyArray_DOUBLE, 1, 0);
if (data == NULL) {
return NULL;
}
npts = data->dimensions[data->nd-1];
data->dimensions[data->nd - 1] = npts/2 + 1;
ret = (PyArrayObject *)PyArray_Zeros(data->nd, data->dimensions,
PyArray_DescrFromType(PyArray_CDOUBLE), 0);
data->dimensions[data->nd - 1] = npts;
rstep = (ret->dimensions[ret->nd - 1])*2;
descr = PyArray_DescrFromType(PyArray_DOUBLE);
if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
goto fail;
}
if (data == NULL || ret == NULL) {
goto fail;
}
if (nsave != npts*2+15) {
PyErr_SetString(ErrorObject, "invalid work array for fft size");
goto fail;
}
nrepeats = PyArray_SIZE(data)/npts;
NPY_SIGINT_ON;
Py_BEGIN_ALLOW_THREADS
#pragma omp parallel for private(i,rptr,dptr) default(shared)
for (i = 0; i < nrepeats; i++) {
rptr = (double *)ret->data + i*rstep;
dptr = (double *)data->data + i*npts;
memcpy((char *)(rptr+1), dptr, npts*sizeof(double));
rfftf(npts, rptr+1, wsave);
rptr[0] = rptr[1];
rptr[1] = 0.0;
}
Py_END_ALLOW_THREADS
NPY_SIGINT_OFF;
PyArray_Free(op2, (char *)wsave);
Py_DECREF(data);
return (PyObject *)ret;
fail:
PyArray_Free(op2, (char *)wsave);
Py_XDECREF(data);
Py_XDECREF(ret);
return NULL;
}
static char fftpack_rfftb__doc__[] ="";
PyObject *
fftpack_rfftb(PyObject *NPY_UNUSED(self), PyObject *args)
{
PyObject *op1, *op2;
PyArrayObject *data, *ret;
PyArray_Descr *descr;
double *wsave, *dptr, *rptr;
npy_intp nsave;
int npts, nrepeats, i;
if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
return NULL;
}
data = (PyArrayObject *)PyArray_ContiguousFromObject(op1,
PyArray_CDOUBLE, 1, 0);
if (data == NULL) {
return NULL;
}
npts = data->dimensions[data->nd - 1];
ret = (PyArrayObject *)PyArray_Zeros(data->nd, data->dimensions,
PyArray_DescrFromType(PyArray_DOUBLE), 0);
descr = PyArray_DescrFromType(PyArray_DOUBLE);
if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
goto fail;
}
if (data == NULL || ret == NULL) {
goto fail;
}
if (nsave != npts*2 + 15) {
PyErr_SetString(ErrorObject, "invalid work array for fft size");
goto fail;
}
nrepeats = PyArray_SIZE(ret)/npts;
NPY_SIGINT_ON;
Py_BEGIN_ALLOW_THREADS
#pragma omp parallel for private(i,rptr,dptr) default(shared)
for (i = 0; i < nrepeats; i++) {
rptr = (double *)ret->data + i*npts;
dptr = (double *)data->data + 2*i*npts;
memcpy((char *)(rptr + 1), (dptr + 2), (npts - 1)*sizeof(double));
rptr[0] = dptr[0];
rfftb(npts, rptr, wsave);
}
Py_END_ALLOW_THREADS
NPY_SIGINT_OFF;
PyArray_Free(op2, (char *)wsave);
Py_DECREF(data);
return (PyObject *)ret;
fail:
PyArray_Free(op2, (char *)wsave);
Py_XDECREF(data);
Py_XDECREF(ret);
return NULL;
}
static char fftpack_rffti__doc__[] ="";
static PyObject *
fftpack_rffti(PyObject *NPY_UNUSED(self), PyObject *args)
{
PyArrayObject *op;
npy_intp dim;
long n;
if (!PyArg_ParseTuple(args, "l", &n)) {
return NULL;
}
/*Magic size needed by rffti*/
dim = 2*n + 15;
/*Create a 1 dimensional array of dimensions of type double*/
op = (PyArrayObject *)PyArray_SimpleNew(1, &dim, PyArray_DOUBLE);
if (op == NULL) {
return NULL;
}
NPY_SIGINT_ON;
Py_BEGIN_ALLOW_THREADS
rffti(n, (double *)((PyArrayObject*)op)->data);
Py_END_ALLOW_THREADS
NPY_SIGINT_OFF;
return (PyObject *)op;
}
/* List of methods defined in the module */
static struct PyMethodDef fftpack_methods[] = {
{"cfftf", fftpack_cfftf, 1, fftpack_cfftf__doc__},
{"cfftb", fftpack_cfftb, 1, fftpack_cfftb__doc__},
{"cffti", fftpack_cffti, 1, fftpack_cffti__doc__},
{"rfftf", fftpack_rfftf, 1, fftpack_rfftf__doc__},
{"rfftb", fftpack_rfftb, 1, fftpack_rfftb__doc__},
{"rffti", fftpack_rffti, 1, fftpack_rffti__doc__},
{NULL, NULL} /* sentinel */
};
/* Initialization function for the module (*must* be called initfftpack) */
static char fftpack_module_documentation[] =
""
;
PyMODINIT_FUNC initfftpack_lite(void)
{
PyObject *m, *d;
/* Create the module and add the functions */
m = Py_InitModule4("fftpack_lite", fftpack_methods,
fftpack_module_documentation,
(PyObject*)NULL,PYTHON_API_VERSION);
/* Import the array object */
import_array();
/* Add some symbolic constants to the module */
d = PyModule_GetDict(m);
ErrorObject = PyErr_NewException("fftpack.error", NULL, NULL);
PyDict_SetItemString(d, "error", ErrorObject);
/* XXXX Add constants here */
}_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
