Module shallow_water::stafft

Fourier transform module. This is not a general purpose transform package but is designed to be quick for arrays of length 2^n. It will work if the array length is of the form 2^i * 3^j * 4^k * 5^l * 6^m (integer powers obviously).

Minimal error-checking is performed by the code below. The only check is that the initial factorisation can be performed. Therefore if the transforms are called with an array of length <2, or a trig array not matching the length of the array to be transformed the code will fail in a spectacular way (eg. Seg. fault or nonsense returned). It is up to the calling code to ensure everything is called sensibly. The reason for stripping error checking is to speed up the backend by performing less if() evaluations - as errors in practice seem to occur very rarely. So the good news is this should be a fast library - the bad is that you may have to pick around in it if there are failures.

To initialise the routines call init(n,factors,trig,ierr). This fills a factorisation array (factors), and a sin/cos array (trig). These must be kept in memory by the calling program. The init routine can be called multiple times with different arrays if more than one length of array is to be transformed. If a factorisation of the array length n cannot be found (as specified above) then the init routine will exit immediately and the integer ierr will be set to 1. If the init returns with ierr=0 then the call was successful.

Top-level subroutines contained in this module are:

1. initfft(n,factors,trig) : Performs intialisation of the module, by working out the factors of n (the FFT length). This will fail if n is not factorised completely by 2,3,4,5,6. The trig array contains the necessary cosine and sine values. Both arrays passed to init must be kept between calls to routines in this module.
2. forfft(m,n,x,trig,factors) : This performs a FFT of an array x containing m vectors of length n. The transform length is n. This inverse of this transform is obtained by revfft.
3. revfft(m,n,x,trig,factors) : This performs an inverse FFT of an array x containing m vectors of length n. The transform length is n. This inverse of this transform is forfft.
4. dct(m,n,x,trig,factors) : This performs a discrete cosine transform of an array x containing m vectors of length n. The transform length is n. This routine calls forfft and performs pre- and post- processing to obtain the transform. This transform is it's own inverse.
5. dst(m,n,x,trig,factors) : This performs a discrete sine transform of an array x containing m vectors of length n. The transform length is n. This routine calls forfft and performs pre- and post- processing to obtain the transform. This transform is it's own inverse.

The storage of the transformed array is in 'Hermitian form'. This means that, for the jth vector the values x(j,1:nw) contain the cosine modes of the transform, while the values x(j,nw+1:n) contain the sine modes (in reverse order ie. wave number increasing from n back to nw+1). [Here, for even n, nw=n/2, and for odd n, nw=(n-1)/2].

Functions

factorisen
forfft	Main physical to spectral (forward) FFT routine. Performs m transforms of length n in the array x which is dimensioned x(m,n). The arrays trig and factors are filled by the init routine and should be kept from call to call. Backend consists of mixed-radix routines, with 'decimation in time'. Transform is stored in Hermitian form.
forrdx2	Radix two physical to Hermitian FFT with 'decimation in time'.
forrdx3	Radix three physical to Hermitian FFT with 'decimation in time'.
forrdx4	Radix four physical to Hermitian FFT with 'decimation in time'.
forrdx5	Radix five physical to Hermitian FFT with 'decimation in time'.
forrdx6	Radix six physical to Hermitian FFT with 'decimation in time'.
initfft	Subroutine performs initialisation work for all the transforms. It calls routines to factorise the array length n and then sets up a trig array full of sin/cos values used in the transform backend.
revfft	Main spectral to physical (reverse) FFT routine. Performs m reverse transforms of length n in the array x which is dimensioned x(m,n). The arrays trig and factors are filled by the init routine and should be kept from call to call. Backend consists of mixed-radix routines, with 'decimation in frequency'. Reverse transform starts in Hermitian form.
revrdx2	Radix two Hermitian to physical FFT with 'decimation in frequency'.
revrdx3	Radix three Hermitian to physical FFT with 'decimation in frequency'.
revrdx4	Radix four Hermitian to physical FFT with 'decimation in frequency'.
revrdx5	Radix five Hermitian to physical FFT with 'decimation in frequency'.
revrdx6	Radix six Hermitian to physical FFT with 'decimation in frequency'.