[PYTHON MATRIX-SIG] fft's

[David Ascher]

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
=================