[PYTHON MATRIX-SIG] fft's
David Ascher
da@maigret.cog.brown.edu
Wed, 28 Feb 1996 10:54:59 -0500 (EST)
I've been working on python interface to the FFT routines in the
complib.sgimath library (since those are the ones I have optimized
binaries for).
This has been educational, but it has raised several questions which I
can't answer. First I'll talk about a specific one related to 2-d and
3-d FFT's, then ask a few questions about interfacing to fortran in
general.
The complib.sgimath (based on stuff from netlib) has a set of fft
calls, 1, 2 and 3-d. Things are working fine in 1-D, but I have to
admit that my grokking of 2-D and 3-D FFT's is rather limited.
When doing a 1-d fft, the output array is bigger than the input array by
enough to fit an aligned complex. Fine.
The 2-d fft calls in the fft library I'm using have the following syntax:
int dfft2du( int job, int n1, int n2, double *sequence,
int lda, double *workspace)
job specifies direction (forward vs. backward)
sequence is where the data is stored (both before and after the fft)
workspace is the sines & cosines and the factorization of n1 and n2.
All that is fine.
n1 and n2 are defined as follows:
N1 - INTEGER.
On entry, N1 specifies the number of elements in the
first dimension of the sequence.
Unchanged on exit.
N2 - INTEGER.
On entry, N2 specifies the number of elements in the
second dimension of the sequence.
Unchanged on exit.
and LDA is defined as:
LDA - INTEGER.
On entry, LDA specifies the first dimension of the
array SEQUENCE as declared in the calling (sub)program.
LDA must be at least max( 1, 2*((N1+2)/2) ).
Unchanged on exit.
This is where I get lost. I assumed that if I had a n1 x n2 array to
start with, I'd have an array which would be slightly bigger in, say, n1
and n2 (enough to store a complex?). But this LDA parameter is
completely confusing me. I have similar problems in 3d of course, where
LDA is replaced by ld1 and ld2. I suspect if I grok LDA i'll grok LD1
and LD2. =)
I realize that fortran and C/Python index into arrays differently (row
vs. column first). What is a good way of dealing with this?
This brings to mind a different issue regarding interfaces to existing
libraries. The library I'm using returns rank 1 arrays, which are
really rank 1, 2, or 3 but designed for antiquated languages. Should I
do the conversion to the appropriate rank (e.g. the assignment of shape)
in a .py wrapper, so that the library module itself is as compatible w/
the fortran and C interfaces, or should I make the library
python-friendly, and have it return appropriately sized arrays?
Thanks for reading this far.
--david
=================
MATRIX-SIG - SIG on Matrix Math for Python
send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
=================