CCFFT3D(3S)CCFFT3D(3S)NAME
CCFFT3D, ZZFFT3D - Applies a three-dimensional complex-to-complex Fast
Fourier Transform (FFT)
SYNOPSIS
Single precision complex -> Single precision complex
Fortran:
CALL CCFFT3D (isign, n1, n2, n3, scale, x, ldx, ldx2, y, ldy,
ldy2, table, work, isys)
C/C++:
#include <scsl_fft.h>
int ccfft3d (int isign, int n1, int n2, int n3, float scale,
scsl_complex *x, int ldx, int ldx2, scsl_complex *y, int ldy,
int ldy2, float *table, float *work, int *isys);
C++ STL:
#include <complex.h>
#include <scsl_fft.h>
int ccfft3d (int isign, int n1, int n2, int n3, float scale,
complex<float> *x, int ldx, int ldx2, complex<float> *y, int
ldy, int ldy2, float *table, float *work, int *isys);
Double precision complex -> Double precision complex
Fortran:
CALL ZZFFT3D (isign, n1, n2, n3, scale, x, ldx, ldx2, y, ldy,
ldy2, table, work, isys)
C/C++:
#include <scsl_fft.h>
int zzfft3d (int isign, int n1, int n2, int n3, double scale,
scsl_zomplex *x, int ldx, int ldx2, scsl_zomplex *y, int ldy,
int ldy2, double *table, double *work, int *isys);
C++ STL:
#include <complex.h>
#include <scsl_fft.h>
int zzfft3d (int isign, int n1, int n2, int n3, double scale,
complex<double> *x, int ldx, int ldx2, complex<double> *y, int
ldy, int ldy2, double *table, double *work, int *isys);
IMPLEMENTATION
These routines are part of the SCSL Scientific Library and can be loaded
using either the -lscs or the -lscs_mp option. The -lscs_mp option
directs the linker to use the multi-processor version of the library.
When linking to SCSL with -lscs or -lscs_mp, the default integer size is
4 bytes (32 bits). Another version of SCSL is available in which integers
are 8 bytes (64 bits). This version allows the user access to larger
memory sizes and helps when porting legacy Cray codes. It can be loaded
Page 1
CCFFT3D(3S)CCFFT3D(3S)
by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
only one of the two versions; 4-byte integer and 8-byte integer library
calls cannot be mixed.
The C and C++ prototypes shown above are appropriate for the 4-byte
integer version of SCSL. When using the 8-byte integer version, the
variables of type int become long long and the <scsl_fft_i8.h> header
file should be included.
DESCRIPTION
These routines compute the three-dimensional complex FFT of the complex
matrix X, and store the results in the complex matrix Y.
In FFT applications, it is customary to use zero-based subscripts; the
formulas are simpler that way. So suppose the arrays are declared as
follows:
Fortran:
COMPLEX X(0:ldx-1, 0:ldx2-1, 0:n3-1)
COMPLEX Y(0:ldy-1, 0:ldy2-1, 0:n3-1)
C/C++:
scsl_complex x[n3][ldx2][ldx], y[n3][ldy2][ldy];
C++ STL:
complex<float> x[n3][ldx2][ldx], y[n3][ldy2][ldy];
These routines compute the formula:
Y(k1,k2,k3) =
n1-1 n2-1 n3-1
scale*Sum Sum Sum [X(j1,j2,j3)*w1**(j1*k1)*w2**(j2*k2)*w3**(j3*k3)]
j1=0 j2=0 j3=0
for k1 = 0, ..., n1 - 1,
k2 = 0, ..., n2 - 1,
k3 = 0, ..., n3 - 1,
where:
w1 = exp(isign*2*pi*i/n1),
Page 2
CCFFT3D(3S)CCFFT3D(3S)
w2 = exp(isign*2*pi*i/n2),
w3 = exp(isign*2*pi*i/n3),
i = + sqrt(-1)
pi = 3.14159...
isign = +1 or -1
Different authors use different conventions for which of the transforms,
isign = +1 or isign = -1, is the forward or inverse transform, and what
the scale factor should be in either case. You can make this routine
compute any of the various possible definitions, however, by choosing the
appropriate values for isign and scale.
The relevant fact from FFT theory is this: If you take the FFT with any
particular values of isign and scale, the mathematical inverse function
is computed by taking the FFT with -isign and 1/(n1 * n2 * n3 * scale).
In particular, if you use isign = +1 and scale = 1.0 for the forward FFT,
you can compute the inverse FFT by using isign = -1 and scale = 1/(n1 *
n2 * n3).
See the NOTES section of this man page for information about the
interpretation of the data types described in the following arguments.
This routine has the following arguments:
isign Integer. (input)
Specifies whether to initialize the table array or to do the
forward or inverse Fourier transform, as follows:
If isign = 0, the routine initializes the table array and
returns. In this case, the only arguments used or checked are
isign, n1, n2, n3, and table.
If isign = +1 or -1, the value of isign is the sign of the
exponent used in the FFT formula.
n1 Integer. (input)
Transform size in the first dimension. If n1 is not positive,
the routine returns without computing a transform.
n2 Integer. (input)
Transform size in the second dimension. If n2 is not positive,
the routine returns without computing a transform.
n3 Integer. (input)
Transform size in the third dimension. If n3 is not positive,
the routine returns without computing a transform.
Page 3
CCFFT3D(3S)CCFFT3D(3S)
scale Scale factor. (input)
CCFFT3D: Single precision.
ZZFFT3D: Double precision.
Each element of the output array is multiplied by scale after
taking the Fourier transform, as defined previously.
x Array of dimensions (ldx, ldx2, n3). (input)
CCFFT3D: Single precision complex array.
ZZFFT3D: Double complex array.
Input array of values to be transformed.
ldx Integer. (input)
The first dimension of x, as it was declared in the calling
program (the leading dimension of x). ldx >= MAX(n1, 1).
ldx2 Integer. (input)
The second dimension of x, as it was declared in the calling
program. ldx2 >= MAX(n2, 1).
y Array of dimensions (ldy, ldy2, n3). (output)
CCFFT3D: Single precision complex array.
ZZFFT3D: Double complex array.
Output array of transformed values. The output array may be
the same as the input array, in which case, the transform is
done in place; that is, the input array is overwritten with the
transformed values. In this case, it is necessary that
ldx = ldy, and ldx2 = ldy2.
ldy Integer. (input)
The first dimension of y, as it was declared in the calling
program (the leading dimension of y). ldy >= MAX(n1, 1).
ldy2 Integer. (input)
The second dimension of y, as it was declared in the calling
program. ldy2 >= MAX(n2, 1).
table Array of dimension (2* n1 + NF) + (2 * n2 + NF) + (2 * n3 + NF)
(input or output)
CCFFT3D: Single precision array.
ZZFFT3D: Double precision array.
Table of factors and root of unity. See the description of the
isys argument for the value of NF.
If isign = 0, the routine initializes table (table is output
only).
If isign = +1 or -1, the values in table are assumed to be
initialized already by a prior call with isign = 0 (table is
input only).
Page 4
CCFFT3D(3S)CCFFT3D(3S)
work Array of dimension 2 * MAX(n1, n2, n3)
CCFFT3D: Single precision array.
ZZFFT3D: Double precision array.
Work array. This is a scratch array used for intermediate
calculations. Its address space must be different from that of
the input and output arrays.
isys Integer array dimensioned 0..isys(0).
An array that gives implementation-specific information. All
features and functions of the FFT routines specific to any
particular implementation are confined to this isys array.
In the Origin series implementation, isys(0)=0 and isys(0)=1
are supported. In SCSL versions prior to 1.3, only isys(0)=0
was allowed. For isys(0)=0, NF=30, and for isys(0)=1, NF=256.
The NF words of storage in the table array contain a
factorization of the length of the transform.
The smaller value of NF for isys(0)=0 is historical. It is too
small to store all the required factors for the highest
performing FFT, so when isys(0)=0, extra space is allocated
when the table array is initialized. To avoid memory leaks,
this extra space must be deallocated when the table array is no
longer needed. The CCFFT3DF routine is used to release this
memory. Due to the potential for memory leaks, the use of
isys(0)=0 should be avoided.
For isys(0)=1, the value of NF is large enough so that no extra
memory needs to be allocated, and there is no need to call
CCFFT3DF to release memory. If called, it does nothing.
NOTE: isys(0)=1 means that isys is an integer array with two
elements. The second element, isys(1), will not be accessed.
NOTES
The following data types are described in this documentation:
Term Used Data type
Fortran:
Array dimensioned 0..n-1 x(0:n-1)
Array of dimensions (m,n) x(m,n)
Array of dimensions (m,n,p) x(m,n,p)
Integer INTEGER (INTEGER*8 for -lscs_i8[_mp])
Page 5
CCFFT3D(3S)CCFFT3D(3S)
Single precision REAL
Double precision DOUBLE PRECISION
Single precision complex COMPLEX
Double precision complex DOUBLE COMPLEX
C/C++:
Array dimensioned 0..n-1 x[n]
Array of dimensions (m,n) x[m*n] or x[n][m]
Array of dimensions (m,n,p) x[m*n*p] or x[p][n][m]
Integer int (long long for -lscs_i8[_mp])
Single precision float
Double precision double
Single precision complex scsl_complex
Double precision complex scsl_zomplex
C++ STL:
Array dimensioned 0..n-1 x[n]
Array of dimensions (m,n) x[m*n] or x[n][m]
Array of dimensions (m,n,p) x[m*n*p] or x[p][n][m]
Integer int (long long for -lscs_i8[_mp])
Single precision float
Double precision double
Single precision complex complex<float>
Double precision complex complex<double>
CAUTIONS
Transform sizes with a prime factor exceeding 232-1 are not supported for
the 8-byte integer version of the library.
In addition to the work array, the FFT routines also dynamically allocate
scratch space from the stack. The amount of space allocated can be
slightly bigger than the size of the largest processor cache. For single
processor runs, the default stack size is large enough that these
Page 6
CCFFT3D(3S)CCFFT3D(3S)
allocations generally cause no problems. But for parallel runs, you need
to ensure that the stack size of slave threads is big enough to hold this
scratch space. Failure to reserve sufficient stack space will cause
programs to dump core due to stack overflows. The stack size of MP
library slave threads is controlled via the MP_SLAVE_STACKSIZE
environment variable or the mp_set_slave_stacksize() library routine. See
the mp(3C), mp(3F) and pe_environ(5) reference pages for more information
on controlling the slave stack size. For pthreads applications, the
thread's stack size is specified as one of many creation attributes
provided in the pthread_attr_t argument to pthread_create(3P). The
stacksize attribute should be set explicitly to a non-default value using
the pthread_attr_setstacksize(3P) call, described in the
pthread_attr_init(3P) man page.
Care must be exercised if copies of the table array are used: even though
a copy exists, the original must persist. As an example, the following
code will not work:
#include <scsl_fft.h>
scsl_complex x[129][129][129], y[129][129][129];
float table[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
float work[2*128];
int isys[2];
isys[0] = 1;
{
float table_orig[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
ccfft3d(0, 128, 128, 128, 1.0f, (scsl_complex *) x, 129,
129, (scsl_complex *) y, 129, 129, table_orig,
work, isys);
bcopy(table_orig, table,
((2*128+256)+(2*256+256)+(2*256+256))*sizeof(float));
}
ccfft3d(1, 128, 128, 128, 1.0f, (scsl_complex *) x, 129,
129, (scsl_complex *) y, 129, 129, table, work, isys);
In this example, because table_orig is a stack variable that does not
persist outside of the code block delimited by the braces, the data in
the copy, table, are not guaranteed to be valid. However, the following
code will work because table_orig is persistent:
#include <scsl_fft.h>
scsl_complex x[129][129][129], y[129][129][129];
float table_orig[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
float table[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
float work[2*128];
int isys[2];
isys[0] = 1;
ccfft3d(0, 128, 128, 128, 1.0f, (scsl_complex *) x, 129,
129, (scsl_complex *) y, 129, 129, table_orig,
work, isys);
Page 7
CCFFT3D(3S)CCFFT3D(3S)
bcopy(table_orig, table,
((2*128+256)+(2*256+256)+(2*256+256))*sizeof(float));
ccfft3d(1, 128, 128, 128, 1.0f, (scsl_complex *) x, 129,
129, (scsl_complex *) y, 129, 129, table,
work, isys);
EXAMPLES
The following examples are for Origin series only.
Example 1: Initialize the TABLE array in preparation for doing a three-
dimensional FFT of size 128 by 128 by 128. In this case, only the isign,
n1, n2, n3, and table arguments are used; you can use dummy arguments or
zeros for other arguments.
Fortran:
REAL TABLE ((2*128 + 256) + (2*128 + 256) + (2*128 + 256))
INTEGER ISYS(0:1)
ISYS(0) = 1
CALL CCFFT3D (0, 128, 128, 128, 0.0, DUMMY, 1, 1, DUMMY, 1, 1,
& TABLE, DUMMY, ISYS)
C/C++:
#include <scsl_fft.h>
float table[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
int isys[2];
isys[0] = 1;
ccfft3d(0, 128, 128, 128, 0.0f, NULL, 1, 1, NULL, 1, 1,
table, NULL, isys);
C++ STL:
#include <complex.h>
#include <scsl_fft.h>
float table[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
int isys[2];
isys[0] = 1;
ccfft3d(0, 128, 128, 128, 0.0f, NULL, 1, 1, NULL, 1, 1,
table, NULL, isys);
Example 2: X and Y are complex arrays of dimension (0:128, 0:128,
0:128). The first 128 elements of each dimension contain data; for
performance reasons, the extra element forces the leading dimensions to
be odd numbers. Take the three-dimensional FFT of X and store it in Y.
Initialize the TABLE array, as in example 1.
Page 8
CCFFT3D(3S)CCFFT3D(3S)
Fortran:
COMPLEX X(0:128, 0:128, 0:128)
COMPLEX Y(0:128, 0:128, 0:128)
REAL TABLE((2*128 + 256) + (2*128 + 256) + (2*128 + 256))
REAL WORK (2*128)
INTEGER ISYS(0:1)
ISYS(0) = 1
CALL CCFFT3D(0, 128, 128, 128, 1.0, X, 129, 129, Y, 129, 129,
& TABLE, WORK, ISYS)
CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129, Y, 129, 129,
& TABLE, WORK, ISYS)
C/C++:
#include <scsl_fft.h>
scsl_complex x[129][129][129], y[129][129][129];
float table[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
float work[2*128];
int isys[2];
isys[0] = 1;
ccfft3d(0, 128, 128, 128, 1.0f, (scsl_complex *) x, 129, 129,
(scsl_complex *) y, 129, 129, table, work, isys);
ccfft3d(1, 128, 128, 128, 1.0f, (scsl_complex *) x, 129, 129,
(scsl_complex *) y, 129, 129, table, work, isys);
C++ STL:
#include <complex.h>
#include <scsl_fft.h>
complex<float> x[129][129][129], y[129][129][129];
float table[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
float work[2*128];
int isys[2];
isys[0] = 1;
ccfft3d(0, 128, 128, 128, 1.0f, (complex<float> *) x, 129, 129,
(complex<float> *) y, 129, 129, table, work, isys);
ccfft3d(1, 128, 128, 128, 1.0f, (complex<float> *) x, 129, 129,
(complex<float> *) y, 129, 129, table, work, isys);
Example 3: With X and Y as in example 2, take the inverse FFT of Y and
store it back in X. The scale factor 1.0/(128.0*128.0*128.0) is used.
Assume that the TABLE array is already initialized.
Fortran:
CALL CCFFT3D(-1, 128, 128, 128, 1.0/(128.0**3), Y, 129, 129,
& X, 129, 129, TABLE, WORK, ISYS)
Page 9
CCFFT3D(3S)CCFFT3D(3S)
C/C++:
ccfft3d(-1, 128, 128, 128, 1.0f/(128.0f*128.0f*128.0f),
(scsl_complex *) y, 129, 129, (scsl_complex *) x, 129,
129, table, work, isys);
C++ STL:
ccfft3d(-1, 128, 128, 128, 1.0f/(128.0f*128.0f*128.0f),
(complex<float> *) y, 129, 129, (complex<float> *) x,
129, 129, table, work, isys);
Example 4: Perform the same computation as in example 2, but put the
output back in the array X to save storage space. Use the 8-byte integer
version of SCSL.
Fortran:
COMPLEX X(0:128, 0:128, 0:128)
REAL TABLE((2*128 + 256) + (2*128 + 256) + (2*128 + 256))
REAL WORK (2*128)
INTEGER ISYS(0:1)
ISYS(0) = 1_8
CALL CCFFT3D(0_8, 128_8, 128_8, 128_8, 1.0, X, 129_8, 129_8,
& X, 129_8, 129_8, TABLE, WORK, ISYS)
CALL CCFFT3D(1_8, 128_8, 128_8, 128_8, 1.0, X, 129_8, 129_8,
& X, 129_8, 129_8, TABLE, WORK, ISYS)
C/C++:
#include <scsl_fft_i8.h>
scsl_complex x[129][129][129];
float table[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
float work[2*128];
long long isys[2];
isys[0] = 1LL;
ccfft3d(0LL, 128LL, 128LL, 128LL, 1.0f, (scsl_complex *) x, 129LL,
129LL, (scsl_complex *) x, 129LL, 129LL, table,
work, isys);
ccfft3d(1LL, 128LL, 128LL, 128LL, 1.0f, (scsl_complex *) x, 129LL,
129LL, (scsl_complex *) x, 129LL, 129LL, table,
work, isys);
C++ STL:
#include <complex.h>
#include <scsl_fft_i8.h>
complex<float> x[129][129][129];
Page 10
CCFFT3D(3S)CCFFT3D(3S)
float table[(2*128 + 256) + (2*128 + 256) + (2*128 + 256)];
float work[2*128];
int isys[2];
isys[0] = 1;
ccfft3d(0LL, 128LL, 128LL, 128LL, 1.0f, (complex<float> *) x, 129LL,
129LL, (complex<float> *) x, 129LL, 129LL, table,
work, isys);
ccfft3d(1LL, 128LL, 128LL, 128LL, 1.0f, (complex<float> *) x, 129LL,
129LL, (complex<float> *) x, 129LL, 129LL, table,
work, isys);
Example 5: Perform the same computation as in example 2, but assume that
the lower bound of each Fortran array is 1, rather than 0. The
subroutine calls do not change.
Fortran:
COMPLEX X(129, 129, 129)
COMPLEX Y(129, 129, 129)
...
CALL CCFFT3D(0, 128, 128, 128, 1.0, X, 129, 129,
& Y, 129, 129, TABLE, WORK, ISYS)
CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
& Y, 129, 129, TABLE, WORK, ISYS)
SEE ALSOINTRO_FFT(3S), INTRO_SCSL(3S), CCFFT(3S), CCFFT2D(3S), CCFFTM(3S),
SCFFT(3S), SCFFT2D(3S), SCFFT3D(3S), SCFFTM(3S)
Page 11
