SCSL.h

Classes

SCSL -- C++ Interface to the Sgi/Cray Scientific Library (SCSL) (full description)

class SCSL

Interface

Public Members
static void ccfft(Int isign, Int n, Float scale, Complex* x, Complex* y, Float* table, Float* work, Int isys)
static void ccfft(Int isign, Int n, Double scale, DComplex* x, DComplex* y, Double* table, Double* work, Int isys)
static void zzfft(Int isign, Int n, Double scale, DComplex* x, DComplex* y, Double* table, Double* work, Int isys)
static void scfft(Int isign, Int n, Float scale, Float* x, Complex* y, Float* table, Float* work, Int isys)
static void scfft(Int isign, Int n, Double scale, Double* x, DComplex* y, Double* table, Double* work, Int isys)
static void dzfft(Int isign, Int n, Double scale, Double* x, DComplex* y, Double* table, Double* work, Int isys)
static void csfft(Int isign, Int n, Float scale, Complex* x, Float* y, Float* table, Float* work, Int isys)
static void csfft(Int isign, Int n, Double scale, DComplex* x, Double* y, Double* table, Double* work, Int isys)
static void zdfft(Int isign, Int n, Double scale, DComplex* x, Double* y, Double* table, Double* work, Int isys)
static void ccfftm(Int isign, Int n, Int lot, Float scale, Complex* x, Int ldx, Complex* y, Int ldy, Float* table, Float* work, Int isys)
static void zzfftm(Int isign, Int n, Int lot, Double scale, DComplex* x, Int ldx, DComplex* y, Int ldy, Double* table, Double* work, Int isys)
static void scfftm(Int isign, Int n, Int lot, Float scale, Float* x, Int ldx, Complex* y, Int ldy, Float* table, Float* work, Int isys)
static void dzfftm(Int isign, Int n, Int lot, Double scale, Double* x, Int ldx, DComplex* y, Int ldy, Double* table, Double* work, Int isys)
static void csfftm(Int isign, Int n, Int lot, Float scale, Complex* x, Int ldx, Float* y, Int ldy, Float* table, Float* work, Int isys)
static void zdfftm(Int isign, Int n, Int lot, Double scale, DComplex* x, Int ldx, Double* y, Int ldy, Double* table, Double* work, Int isys)
static void ccfft2d(Int isign, Int n1, Int n2, Float scale, Complex* x, Int ldx, Complex* y, Int ldy, Float* table, Float* work, Int isys)
static void zzfft2d(Int isign, Int n1, Int n2, Double scale, DComplex* x, Int ldx, DComplex* y, Int ldy, Double* table, Double* work, Int isys)
static void scfft2d(Int isign, Int n1, Int n2, Float scale, Float* x, Int ldx, Complex* y, Int ldy, Float* table, Float* work, Int isys)
static void dzfft2d(Int isign, Int n1, Int n2, Double scale, Double* x, Int ldx, DComplex* y, Int ldy, Double* table, Double* work, Int isys)
static void csfft2d(Int isign, Int n1, Int n2, Float scale, Complex* x, Int ldx, Float* y, Int ldy, Float* table, Float* work, Int isys)
static void zdfft2d(Int isign, Int n1, Int n2, Double scale, DComplex* x, Int ldx, Double* y, Int ldy, Double* table, Double* work, Int isys)
static void ccfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Complex* x, Int ldx, Int ldx2, Complex* y, Int ldy, Int ldy2, Float* table, Float* work, Int isys)
static void zzfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, DComplex* x, Int ldx, Int ldx2, DComplex* y, Int ldy, Int ldy2, Double* table, Double* work, Int isys)
static void scfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Float* x, Int ldx, Int ldx2, Complex* y, Int ldy, Int ldy2, Float* table, Float* work, Int isys)
static void dzfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, Double* x, Int ldx, Int ldx2, DComplex* y, Int ldy, Int ldy2, Double* table, Double* work, Int isys)
static void csfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Complex* x, Int ldx, Int ldx2, Float* y, Int ldy, Int ldy2, Float* table, Float* work, Int isys)
static void zdfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, DComplex* x, Int ldx, Int ldx2, Double* y, Int ldy, Int ldy2, Double* table, Double* work, Int isys)

Description

Synopsis

These are C++ wrapper functions for the transform routines in the SGI/Cray Scientific Library (SCSL). The purpose of these definitions is to overload the functions so that C++ users can access the functions in SCSL with identical function names.

Warning Currently, the SCSL is available only on SGI machines.

Member Description

static void ccfft(Int isign, Int n, Float scale, Complex* x, Complex* y, Float* table, Float* work, Int isys)
static void ccfft(Int isign, Int n, Double scale, DComplex* x, DComplex* y, Double* table, Double* work, Int isys)
static void zzfft(Int isign, Int n, Double scale, DComplex* x, DComplex* y, Double* table, Double* work, Int isys)

These routines compute the Fast Fourier Transform (FFT) of the complex vector x, and store the result in vector y. ccfft does the complex-to-complex transform and zzfft does the same for double precision arrays.

In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. Suppose that the arrays are dimensioned as follows:

      COMPLEX X(0:N-1), Y(0:N-1)

The output array is the FFT of the input array, using the following formula for the FFT:

                     n-1
      Y(k) = scale * Sum [ X(j)*w**(isign*j*k) ]    for k = 0, ..., n-1
                     j=0

      where:
      w = exp(2*pi*i/n),
      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/(n*scale). In particular, if you use isign = +1 and scale = 1.0 you can compute the inverse FFT by using isign = -1 and scale = 1.0/n.

The output array may be the same as the input array.

Initialization

The table array stores the trigonometric tables used in calculation of the FFT. You must initialize table by calling the routine with isign = 0 prior to doing the transforms. If the value of the problem size, n, does not change, table does not have to be reinitialized.

Dimensions

In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:

      COMPLEX X(0:N-1)
      COMPLEX Y(0:N-1)

However, if you prefer to use the more customary FORTRAN style with subscripts starting at 1 you do not have to change the calling sequence, as in the following (assuming N > 0):

      COMPLEX X(N)
      COMPLEX Y(N)

Example

These examples use the table and workspace sizes appropriate to the Origin series.

Example 1: Initialize the complex array table in preparation for doing an FFT of size 1024. Only the isign, n, and table arguments are used in this case. You can use dummy arguments or zeros for the other arguments in the subroutine call.

      REAL TABLE(30 + 2048)
      CALL CCFFT(0, 1024, 0.0, DUMMY, DUMMY, TABLE, DUMMY, 0)

Example 2: x and y are complex arrays of dimension (0:1023). Take the FFT of x and store the results in y. Before taking the FFT, initialize the table array, as in example 1.

      COMPLEX X(0:1023), Y(0:1023)
      REAL TABLE(30 + 2048)
      REAL WORK(2048)
      CALL CCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
      CALL CCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)

Example 3: Using the same x and y as in example 2, take the inverse FFT of y and store it back in x. The scale factor 1/1024 is used. Assume that the table array is already initialized.

      CALL CCFFT(-1, 1024, 1.0/1024.0, Y, X, TABLE, WORK, 0)

Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. No change was needed in the subroutine calls.

      COMPLEX X(1024), Y(1024)
      CALL CCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
      CALL CCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)

Example 5: Do the same computation as in example 4, but put the output back in array x to save storage space. Assume that table is already initialized.

      COMPLEX X(1024)
      CALL CCFFT(1, 1024, 1.0, X, X, TABLE, WORK, 0)

Input parameters:

isign
Integer. 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, n, and table.

If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.

n
Integer. Size of the transform (the number of values in the input array). n >= 1.
scale
Scale factor. ccfft: real. zzfft: double precision. Each element of the output array is multiplied by scale after taking the Fourier transform, as defined by the previous formula.
x
Array of dimension (0:n-1). ccfft: complex array. zzfft: double complex array.

Input array of values to be transformed.

isys
Integer. Algorithm used; value dependent on hardware system. Currently, no special options are supported; therefore, you must always specify an isys argument as constant 0.
Output parameters:
y
Array of dimension (0:n-1). ccfft: complex array. zzfft: double complex array. Output array of transformed values. The output array may be the same as the input array. In that case, the transform is done in place and the input array is overwritten with the transformed values.
table
Real array; dimension 2*n+30.

Table of factors and trigonometric functions.

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).

work
Real array; dimension 2*n.

Work array. This is a scratch array used for intermediate calculations. Its address space must be different address space from that of the input and output arrays.

static void scfft(Int isign, Int n, Float scale, Float* x, Complex* y, Float* table, Float* work, Int isys)
static void scfft(Int isign, Int n, Double scale, Double* x, DComplex* y, Double* table, Double* work, Int isys)
static void dzfft(Int isign, Int n, Double scale, Double* x, DComplex* y, Double* table, Double* work, Int isys)
static void csfft(Int isign, Int n, Float scale, Complex* x, Float* y, Float* table, Float* work, Int isys)
static void csfft(Int isign, Int n, Double scale, DComplex* x, Double* y, Double* table, Double* work, Int isys)
static void zdfft(Int isign, Int n, Double scale, DComplex* x, Double* y, Double* table, Double* work, Int isys)

scfft/dzfft computes the FFT of the real array x, and it stores the results in the complex array y. csfft/zdfft computes the corresponding inverse complex-to-real transform.

It is customary in FFT applications to use zero-based subscripts; the formulas are simpler that way. For these routines, suppose that the arrays are dimensioned as follows:

      REAL    X(0:n-1)
      COMPLEX Y(0:n/2)

Then the output array is the FFT of the input array, using the following formula for the FFT:

                     n-1
      Y(k) = scale * Sum [ X(j)*w**(isign*j*k) ]    for k = 0, ..., n/2
                     j=0

      where:
      w = exp(2*pi*i/n),
      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 these routines 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 call scfft with any particular values of isign and scale, the mathematical inverse function is computed by calling csfft with -isign and 1/(n*scale). In particular, if you use isign = +1 and scale = 1.0 in scfft for the forward FFT, you can compute the inverse FFT by using ccfft with isign = -1 and scale = 1.0/n.

Real-to-complex FFTs

Notice in the preceding formula that there are n real input values, and n/2 + 1 complex output values. This property is characteristic of real-to-complex FFTs.

The mathematical definition of the Fourier transform takes a sequence of n complex values and transforms it to another sequence of n complex values. A complex-to-complex FFT routine, such as ccfft, will take n complex input values, and produce n complex output values. In fact, one easy way to compute a real-to-complex FFT is to store the input data in a complex array, then call routine ccfft to compute the FFT. You get the same answer when using the scfft routine.

The reason for having a separate real-to-complex FFT routine is efficiency. Because the input data is real, you can make use of this fact to save almost half of the computational work. The theory of Fourier transforms tells us that for real input data, you have to compute only the first n/2 + 1 complex output values, because the remaining values can be computed from the first half of the values by the simple formula:

      Y(k) = conjg(Y(n-k)) for n/2 <= k <= n-1

where the notation conjgY represents the complex conjugate of y.

In fact, in many applications, the second half of the complex output data is never explicitly computed or stored. Likewise, as explained later, only the first half of the complex data has to be supplied for the complex-to-real FFT.

Another implication of FFT theory is that, for real input data, the first output value, Y(0), will always be a real number; therefore, the imaginary part will always be 0. If n is an even number, Y(n/2) will also be real and thus, have zero imaginary parts.

Complex-to-real FFTs

Consider the complex-to-real case. The effect of the computation is given by the preceding formula, but with X complex and Y real.

Generally, the FFT transforms a complex sequence into a complex sequence. However, in a certain application we may know the output sequence is real. Often, this is the case because the complex input sequence was the transform of a real sequence. In this case, you can save about half of the computational work.

According to the theory of Fourier transforms, for the output sequence, Y, to be a real sequence, the following identity on the input sequence, X, must be true:

      X(k) = conjg(X(n-k)) for n/2 <= k <= n-1

And, in fact, the input values X(k) for k > n/2 need not be supplied; they can be inferred from the first half of the input.

Thus, in the complex-to-real routine, csfft, the arrays can be dimensioned as follows:

      COMPLEX X(0:n/2)
      REAL    Y(0:n-1)

There are n/2 + 1 complex input values and n real output values. Even though only n/2 + 1 input values are supplied, the size of the transform is still n in this case, because implicitly you are using the FFT formula for a sequence of length n.

Another implication of the theory is that X(0) must be a real number (that is, it must have zero imaginary part). Also, if n is even, X(n/2) must also be real. Routine CSFFT assumes that these values are real; if you specify a nonzero imaginary part, it is ignored.

Table Initialization

The table array stores the trigonometric tables used in calculation of the FFT. This table must be initialized by calling the routine with isign = 0 prior to doing the transforms. The table does not have to be reinitialized if the value of the problem size, n, does not change. Because scfft and csfft use the same format for table, either can be used to initialize it (note that CCFFT uses a different table format).

Dimensions

In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared (assuming n > 0):

      REAL    X(0:n-1)
      COMPLEX Y(0:n/2)

No change is needed in the calling sequence; however, if you prefer you can use the more customary Fortran style with subscripts starting at 1, as in the following:

      REAL    X(n)
      COMPLEX Y(n/2 + 1)

Example

These examples use the table and workspace sizes appropriate to Origin series.

Example 1: Initialize the complex array TABLE in preparation for doing an FFT of size 1024. In this case only the arguments isign, n, and table are used. You can use dummy arguments or zeros for the other arguments in the subroutine call.

      REAL TABLE(15 + 1024)
      CALL SCFFT(0, 1024, 0.0, DUMMY, DUMMY, TABLE, DUMMY, 0)

Example 2: X is a real array of dimension (0:1023), and Y is a complex array of dimension (0:512). Take the FFT of X and store the results in Y. Before taking the FFT, initialize the TABLE array, as in example 1.

      REAL X(0:1023)
      COMPLEX Y(0:512)
      REAL TABLE(15 + 1024)
      REAL WORK(1024)
      CALL SCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
      CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)

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/1024 is used. Assume that the TABLE array is initialized already.

      CALL CSFFT(-1, 1024, 1.0/1024.0, Y, X, TABLE, WORK, 0)

Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. The subroutine calls are not changed.

      REAL X(1024)
      COMPLEX Y(513)
      CALL SCFFT(0, 1024, 1.0, X, Y, TABLE, WORK, 0)
      CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)

Example 5: Perform the same computation as in example 4, but equivalence the input and output arrays to save storage space. Assume that the TABLE array is initialized already.

      REAL X(1024)
      COMPLEX Y(513)
      EQUIVALENCE ( X(1), Y(1) )
      CALL SCFFT(1, 1024, 1.0, X, Y, TABLE, WORK, 0)

Input parameters:

isign
Integer. 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, n, and table.

If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.

n
Integer. Size of transform. If n <= 2, scfft/dzfft returns without calculating the transform.
scale
Scale factor. scfft: real. dzfft: double precision. csfft: real. zdfft: double precision. Each element of the output array is multiplied by scale after taking the Fourier transform, as defined by the previous formula.
x
Input array of values to be transformed. scfft: real array of dimension (0:n-1). dzfft: double precision array of dimension (0:n-1). csfft: complex array of dimension (0:n/2). zdfft: double complex array of dimension (0:n/2).
isys
Integer array of dimension (0:isys(0)). Use isys to specify certain processor-specific parameters or options. The first element of the array specifies how many more elements are in the array.

If isys(0) = 0, the default values of such parameters are used. In this case, you can specify the argument value as the scalar integer constant 0. If isys(0) > 0, isys(0) gives the upper bound of the isys array; that is, if il = isys(0), user-specified parameters are expected in isys(1) through isys(il).

Output parameters:
y
Output array of transformed values. scfft: complex array of dimension (0:n/2). dzfft: double complex array of dimension (0:n/2). csfft: real array of dimension (0:n-1). zdfft: double precision array of dimension (0:n-1).

The output array, y, is the FFT of the the input array, x, computed according to the preceding formula. The output array may be equivalenced to the input array in the calling program. Be careful when dimensioning the arrays, in this case, to allow for the fact that the complex array contains two (real) words more than the real array.

table
Real array; dimension n+15.

Table of factors and trigonometric functions.

If isign = 0, the table array is initialized to contain trigonometric tables needed to compute an FFT of size n.

If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0.

work
Real array; dimension n.

Work array used for intermediate calculations. Its address space must be different from that of the input and output arrays.

static void ccfftm(Int isign, Int n, Int lot, Float scale, Complex* x, Int ldx, Complex* y, Int ldy, Float* table, Float* work, Int isys)
static void zzfftm(Int isign, Int n, Int lot, Double scale, DComplex* x, Int ldx, DComplex* y, Int ldy, Double* table, Double* work, Int isys)

ccfftm/zzfftm computes the FFT of each column of the complex matrix x, and stores the results in the columns of complex matrix y.

Suppose the arrays are dimensioned as follows:

         COMPLEX X(0:ldx-1, 0:lot-1)
         COMPLEX Y(0:ldy-1, 0:lot-1)
    
    where ldx >= n, ldy >= n.
    

Then column L of the output array is the FFT of column L of the input array, using the following formula for the FFT:

                        n-1
      Y(k, L) = scale * Sum [ X(j)*w**(isign*j*k) ]
                        j=0
      for k = 0, ..., n-1
          L = 0, ..., lot-1
      where:
          w = exp(2*pi*i/n),
          i = + sqrt(-1),
          pi = 3.14159...,
          isign = +1 or -1
          lot = the number of columns to transform

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/(n * scale). In particular, if you use isign = +1 and scale = 1.0 for the forward FFT, you can compute the inverse FFT by using the following: isign = -1 and scale = 1.0/n.

This section contains information about the algorithm for these routines, the initialization of the table array, the declaration of dimensions for x and y arrays, some performance tips, and some implementation-dependent details.

Algorithm

These routines use decimation-in-frequency type FFT. It takes the FFT of the columns and vectorizes the operations along the rows of the matrix. Thus, the vector length in the calculations depends on the row size, and the strides for vector loads and stores are the leading dimensions, ldx and ldy.

Initialization

The table array stores the trigonometric tables used in calculation of the FFT. You must initialize the table array by calling the routine with isign = 0 prior to doing the transforms. If the value of the problem size, n, does not change, table does not have to be reinitialized.

Dimensions

In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:

      COMPLEX X(0:ldx-1, 0:lot-1)
      COMPLEX Y(0:ldy-1, 0:lot-1)

The calling sequence does not have to change, however, if you prefer to use the more customary Fortran style with subscripts starting at 1. The same values of ldx and ldy would be passed to the subroutine even if the input and output arrays were dimensioned as follows:

      COMPLEX X(ldx, lot)
      COMPLEX Y(ldy, lot)

Example

Example 1: Initialize the TABLE array in preparation for doing an FFT of size 128. Only the isign, n, and table arguments are used in this case. You can use dummy arguments or zeros for the other arguments in the subroutine call.

       REAL TABLE(30 + 256)
       CALL CCFFTM(0, 128, 0, 0., DUMMY, 1, DUMMY, 1, TABLE, DUMMY, 0)

Example 2: X and Y are complex arrays of dimension (0:128) by (0:55). The first 128 elements of each column contain data. For performance reasons, the extra element forces the leading dimension to be an odd number. Take the FFT of the first 50 columns of X and store the results in the first 50 columns of Y. Before taking the FFT, initialize the TABLE array, as in example 1.

       COMPLEX X(0:128, 0:55)
       COMPLEX Y(0:128, 0:55)
       REAL TABLE(30 + 256)
       REAL WORK(256)
       ...
       CALL CCFFTM(0, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
       CALL CCFFTM(1, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)

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/128 is used. Assume that the TABLE array is already initialized.

       CALL CCFFTM(-1, 128, 50, 1./128., Y, 129, X, 129, TABLE,WORK,0)

Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. The subroutine calls are not changed.

       COMPLEX X(129, 55)
       COMPLEX Y(129, 55)
       ...
       CALL CCFFTM(0, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
       CALL CCFFTM(1, 128, 50, 1.0, X, 129, Y, 129, TABLE, WORK, 0)

Example 5: Perform the same computation as in example 4, but put the output back in array X to save storage space. Assume that the TABLE array is already initialized.

       COMPLEX X(129, 55)
       ...
       CALL CCFFTM(1, 128, 50, 1.0, X, 129, X, 129, TABLE, WORK, 0)

Input parameters:

isign
Integer. 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, n, and table.

If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.

n
Integer. Size of each transform (the number of elements in each column of the input and output matrix to be transformed). Performance depends on the value of n, as explained above. n >= 0; if n = 0, the routine returns.
lot
Integer. The number of transforms to be computed (lot size). This is the number of elements in each row of the input and output matrix. lot >= 0. If lot = 0, the routine returns.
scale
Scale factor. ccfftm: real. zzfftm: double precision. Each element of the output array is multiplied by scale factor after taking the Fourier transform, as defined previously.
x
Array of dimension (0:ldx-1, 0:n2-1). ccfftm: real array. zzfftm: double precision array. Input array of values to be transformed.
ldx
The number of rows in the x array, as it was declared in the calling program (the leading dimension of X). ldx >= MAX(n, 1).
ldy
Integer. The number of rows in the y array, as it was declared in the calling program (the leading dimension of y). ldy >= MAX(n, 1).
isys
Integer array of dimension (0:isys(0)). The first element of the array specifies how many more elements are in the array. Use isys to specify certain processor-specific parameters or options.

If isys(0) = 0, the default values of such parameters are used. In this case, you can specify the argument value as the scalar integer constant 0.

If isys(0) > 0, isys(0) gives the upper bound of the isys array. Therefore, if il = isys(0), user-specified parameters are expected in isys(1) through isys(il).

Output parameters:
y
Array of dimension (0:ldy-1, 0:lot-1). ccfftm: complex array. zzfftm: double complex array. Output array of transformed values. Each column of the output array, y, is the FFT of the corresponding column of the input array, x, computed according to the preceding formula.

The output array may be the same as the input array. In that case, the transform is done in place. The input array is overwritten with the transformed values. In this case, it is necessary that ldx = ldy.

table
Real array; dimension (30 + 2n). Table of factors and trigonometric functions.

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).

work
Real array; dimension 2n. 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.

static void scfftm(Int isign, Int n, Int lot, Float scale, Float* x, Int ldx, Complex* y, Int ldy, Float* table, Float* work, Int isys)
static void dzfftm(Int isign, Int n, Int lot, Double scale, Double* x, Int ldx, DComplex* y, Int ldy, Double* table, Double* work, Int isys)
static void csfftm(Int isign, Int n, Int lot, Float scale, Complex* x, Int ldx, Float* y, Int ldy, Float* table, Float* work, Int isys)
static void zdfftm(Int isign, Int n, Int lot, Double scale, DComplex* x, Int ldx, Double* y, Int ldy, Double* table, Double* work, Int isys)

scfftm/dzfftm computes the FFT of each column of the real matrix X, and it stores the results in the corresponding column of the complex matrix Y. csfftm/zdfftm computes the corresponding inverse transforms.

In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. First, the function of scfftm is described. Suppose that the arrays are dimensioned as follows:

         REAL    X(0:ldx-1, 0:lot-1)
         COMPLEX Y(0:ldy-1, 0:lot-1)
    
    where ldx >= n, ldy >= n/2 + 1.
    

Then column L of the output array is the FFT of column L of the input array, using the following formula for the FFT:

                       n-1
    Y(k, L) = scale *  Sum  [ X(j, L)*w**(isign*j*k) ]
                       j=0
    
    for k = 0, ..., n/2
        L = 0, ..., lot-1 where:
        w = exp(2*pi*i/n),
        i = + sqrt(-1)
        pi = 3.14159...,
        isign = +1 or -1,
        lot = the number of columns to transform
    

Different authors use different conventions for which transform (isign = +1 or isign = -1) is used in the real-to-complex case, and what the scale factor should be. Some adopt the convention that isign = 1 for the real-to-complex transform, and isign = -1 for the complex-to-real inverse. Others use the opposite convention. You can make these routines 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 use scfftm to take the real-to-complex FFT, using any particular values of isign and scale, the mathematical inverse function is computed by using csfftm with -isign and 1/ (n*scale). In particular, if you call scfftm with isign = +1 and scale = 1.0, you can use csfftm to compute the inverse complex-to-real FFT by using isign = -1 and scale = 1.0/n.

Real-to-complex FFTs

Notice in the preceding formula that there are n real input values and (n/2) + 1 complex output values for each column. This property is characteristic of real-to-complex FFTs.

The mathematical definition of the Fourier transform takes a sequence of n complex values and transforms it to another sequence of n complex values. A complex-to-complex FFT routine, such as ccfftm, will take n complex input values and produce n complex output values. In fact, one easy way to compute a real-to-complex FFT is to store the input data x in a complex array, then call routine ccfftm to compute the FFT. You get the same answer when using the scfftm routine.

A separate real-to-complex FFT routine is more efficient than the equivalent complex-to-complex routine. Because the input data is real, you can make use of this fact to save almost half of the computational work. According to the theory of Fourier transforms, for real input data, you have to compute only the first n/2 + 1 complex output values in each column, because the second half of the FFT values in each column can be computed from the first half of the values by the simple formula:

         Y    = conjgY          for n/2 <= k <= n-1
          k,L          n-k, L
    
    where the notation conjg(z) represents the complex conjugate of z.
    

In fact, in many applications, the second half of the complex output data is never explicitly computed or stored. Likewise, you must supply only the first half of the complex data in each column has to be supplied for the complex-to-real FFT.

Another implication of FFT theory is that for real input data, the first output value in each column, Y(0, L), will always be a real number; therefore, the imaginary part will always be 0. If n is an even number, Y(n/2, L) will also be real and have 0 imaginary parts.

Complex-to-real FFTs

Consider the complex-to-real case. The effect of the computation is given by the preceding formula, but with X complex and Y real.

In general, the FFT transforms a complex sequence into a complex sequence; however, in a certain application you may know the output sequence is real, perhaps because the complex input sequence was the transform of a real sequence. In this case, you can save about half of the computational work.

According to the theory of Fourier transforms, for the output sequence, Y, to be a real sequence, the following identity on the input sequence, X, must be true:

         X    = conjgX         for n/2 <= k <= n-1
          k,L          n-k,L
    And, in fact, the following input values
    
         X    for k > n/2
          k,L
    do not have to be supplied, because they can be inferred from the
    first half of the input.
    

Thus, in the complex-to-real routine, CSFFTM, the arrays can be dimensioned as follows:

         COMPLEX X(0:ldx-1, 0:lot-1)
         REAL    Y(0:ldy-1, 0:lot-1)
    
    where ldx >= n/2 + 1, ldy >= n.
    

In each column, there are (n/2) + 1 complex input values and n real output values. Even though only (n/2) + 1 input values are supplied, the size of the transform is still n in this case, because implicitly the FFT formula for a sequence of length n is used.

Another implication of the theory is that X(0, L) must be a real number (that is, must have zero imaginary part). If n is an even number, X(n/2, L) must also be real. Routine CSFFTM assumes that each of these values is real; if a nonzero imaginary part is given, it is ignored.

Table Initialization

The table array contains the trigonometric tables used in calculation of the FFT. You must initialize this table by calling the routine with isign = 0 prior to doing the transforms. table does not have to be reinitialized if the value of the problem size, n, does not change.

Dimensions

In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared (for SCFFTM):

      REAL    X(0:ldx-1, 0:lot-1)
      COMPLEX Y(0:ldy-1, 0:lot-1)

No change is made in the calling sequence, however, if you prefer to use the more customary Fortran style with subscripts starting at 1. The same values of ldx and ldy would be passed to the subroutine even if the input and output arrays were dimensioned as follows:

      REAL    X(ldx, lot)
      COMPLEX Y(ldy, lot)

Example 1: Initialize the complex array TABLE in preparation for doing an FFT of size 128. In this case only the isign, n, and table arguments are used; you may use dummy arguments or zeros for the other arguments in the subroutine call.

       REAL TABLE(15 + 128)
       CALL SCFFTM(0, 128, 1, 0.0, DUMMY, 1, DUMMY, 1,
      &  TABLE, DUMMY, 0)

Example 2: X is a real array of dimension (0:128, 0:55), and Y is a complex array of dimension (0:64, 0:55). The first 128 elements in each column of X contain data; the extra element forces an odd leading dimension. Take the FFT of the first 50 columns of X and store the results in the first 50 columns of Y. Before taking the FFT, initialize the TABLE array, as in example 1.

       REAL    X(0:128, 0:55)
       COMPLEX Y(0:64,  0:55)
       REAL    TABLE(15 + 128)
       REAL    WORK((128)
       ...
       CALL SCFFTM(0, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
       CALL SCFFTM(1, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)

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/128 is used. Assume that the TABLE array is initialized already.

       CALL CSFFTM(-1, 128, 50, 1.0/128.0, Y, 65, X, 129,
      &  TABLE, WORK, 0)

Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. No change is made in the subroutine calls.

       REAL    X(129, 56)
       COMPLEX Y(65, 56)
       ...
       CALL SCFFTM(0, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
       CALL SCFFTM(1, 128, 50, 1.0, X, 129, Y, 65, TABLE, WORK, 0)

Example 5: Perform the same computation as in example 4, but equivalence the input and output arrays to save storage space. In this case, a row must be added to X, because it is equivalenced to a complex array. The leading dimension of X is two times an odd number; therefore, memory bank conflicts are minimal. Assume that TABLE is initialized already.

       REAL    X(130, 56)
       COMPLEX Y(65, 56)
       EQUIVALENCE ( X(1, 1), Y(1, 1) )
       ...
       CALL SCFFTM(1, 128, 50, 1.0, X, 130, Y, 65, TABLE, WORK, 0)

Input parameters:

isign
Integer. 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, n, and table.

If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.

n
Integer. Size of the transforms (the number of elements in each column of the input and output matrix to be transformed). If n is not positive, scfftm or csfftm returns without computing a transforms.
lot
Integer. The number of transforms to be computed (or "lot size"). This is the number of elements in each row of the input and output matrix. If lot is not positive, csfftm or scfftm returns without computing a transforms.
scale
Scale factor. scfftm: real. dzfftm: double precision. csfftm: real. zdfftm: double precision. Each element of the output array is multiplied by scale after taking the transform, as defined in the preceding formula.
x
Input array of values to be transformed. Dimension (0:ldx-1, 0:lot-1). scfftm: real array. dzfftm: double precision array. csfftm: complex array. zdfftm: double complex array.
ldx
Integer. The number of rows in the x array, as it was declared in the calling program. That is, the leading dimension of x. scfftm, dzfftm: ldx >= MAX(n, 1). csfftm, zdfftm: ldx >= MAX(n/2 + 1, 1).
ldy
Integer. The number of rows in the y array, as it was declared in the calling program (the leading dimension of y). scfftm, dzfftm: ldy >= MAX(n/2 + 1, 1). csfftm, zdfftm: ldy >= MAX(n, 1).
isys
Integer array of dimension (0:isys(0)). The first element of the array specifies how many more elements are in the array. Use isys to specify certain processor-specific parameters or options.

If isys(0) = 0, the default values of such parameters are used. In this case, you can specify the argument value as the scalar integer constant 0.

If isys(0) > 0, isys(0) gives the upper bound of the isys array. Therefore, if il = isys(0), user-specified parameters are expected in isys(1) through isys(il).

Output parameters:
y
Output array of transformed values. Dimension (0:ldy-1, 0:lot-1). scfftm: complex array. dzfftm: double complex array. csfftm: real array. zdfftm: double precision array.

Each column of the output array, y, is the FFT of the corresponding column of the input array, x, computed according to the preceding formula. The output array may be equivalenced to the input array. In that case, the transform is done in place and the input array is overwritten with the transformed values. In this case, the following conditions on the leading dimensions must hold:

scfftm, dzfftm: ldx = 2ldy. csfftm, zdfftm: ldy = 2ldx.

table
Real array; dimension (15 + n). Table of factors and trigonometric functions. This array must be initialized by a call to scfftm (or csfftm) with isign = 0.

If isign = 0, table is initialized to contain trigonometric tables needed to compute an FFT of length n.

work
Real array; dimension n. Work array used for intermediate calculations. Its address space must be different from that of the input and output arrays.

static void ccfft2d(Int isign, Int n1, Int n2, Float scale, Complex* x, Int ldx, Complex* y, Int ldy, Float* table, Float* work, Int isys)
static void zzfft2d(Int isign, Int n1, Int n2, Double scale, DComplex* x, Int ldx, DComplex* y, Int ldy, Double* table, Double* work, Int isys)

These routines compute the two-dimensional complex Fast Fourier Transform (FFT) of the complex matrix x, and store the results in the complex matrix y. ccfft2d does the complex-to-complex transform and zzfft does the same for double precision arrays.

In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. Suppose that the arrays are dimensioned as follows:

      COMPLEX X(0:n1-1, 0:n2-1)
      COMPLEX Y(0:n1-1, 0:n2-1)

These routines compute the formula:

                        n2-1  n1-1
    Y(k1, k2) = scale * Sum   Sum [ X(j1, j2)*w1**(j1*k1)*w2**(j2*k2) ]
                        j2=0  j1=0

    for k1 = 0, ..., n1-1
        k2 = 0, ..., n2-1

    where:
        w1 = exp(isign*2*pi*i/n1)
        w2 = exp(isign*2*pi*i/n2)
        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*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.0/(n1*n2).

Algorithm

These routines use a routine very much like ccfftm/zzfftm to do multiple FFTs first on all columns in an input matrix and then on all of the rows.

Initialization

The table array stores factors of n1 and n2 and also trigonometric tables that are used in calculation of the FFT. This table must be initialized by calling the routine with isign = 0. If the values of the problem sizes, n1 and n2, do not change, the table does not have to be reinitialized.

Dimensions

In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:

      COMPLEX X(0:ldx-1, 0:n2-1)
      COMPLEX Y(0:ldy-1, 0:n2-1)

However, the calling sequence does not change if you prefer to use the more customary Fortran style with subscripts starting at 1. The same values of ldx and ldy would be passed to the subroutine even if the input and output arrays were dimensioned as follows:

      COMPLEX X(ldx, n2)
      COMPLEX Y(ldy, n2)

Example

All examples here are for Origin series only.

Example 1: Initialize the TABLE array in preparation for doing a two-dimensional FFT of size 128 by 256. In this case only the isign, n1, n2, and table arguments are used; you can use dummy arguments or zeros for other arguments.

        REAL TABLE ((30 + 256) + (30 + 512))
        CALL CCFFT2D (0, 128, 256, 0.0, DUMMY, 1, DUMMY, 1,
       &  TABLE, DUMMY, 0)

Example 2: X and Y are complex arrays of dimension (0:128, 0:255). The first 128 elements of each column contain data. For performance reasons, the extra element forces the leading dimension to be an odd number. Take the two-dimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.

       COMPLEX X(0:128, 0:255)
       COMPLEX Y(0:128, 0:255)
       REAL     TABLE((30 + 256) + (30 + 512))
       REAL     WORK(2*128*256)
       ...
       CALL CCFFT2D(0, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
       CALL CCFFT2D(1, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)

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/(128*256) is used. Assume that the TABLE array is already initialized.

       CALL CCFFT2D(-1, 128, 256, 1.0/(128.0*256.0), Y, 129,
      &  X, 129, TABLE, WORK, 0)

Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. The subroutine calls are not changed.

       COMPLEX X(129, 256)
       COMPLEX Y(129, 256)
       ...
       CALL CCFFT2D(0, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)
       CALL CCFFT2D(1, 128, 256, 1.0, X, 129, Y, 129, TABLE, WORK, 0)

Example 5: Perform the same computation as in example 4, but put the output back in array X to save storage space. Assume that the TABLE array is already initialized.

       COMPLEX X(129, 256)
       ...
       CALL CCFFT2D(1, 128, 256, 1.0, X, 129, X, 129, TABLE, WORK, 0)

Input parameters:

isign
Integer. Specifies whether to initialize the table array or to do the forward or inverse 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, table.

If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.

n1
Integer. Transform size in the first dimension. If n1 is not positive, the routine returns without performing a transform.
n2
Integer. Transform size in the second dimension. If n2 is not positive, the routine returns without performing a transform.
scale
Scale factor. ccfft2d: real. zzfft2d: double precision. Each element of the output array is multiplied by scale factor after taking the Fourier transform, as defined previously.
x
Array of dimension (0:ldx-1, 0:n2-1). ccfft2d: complex array. zzfft2d: double complex array. Input array of values to be transformed.
ldx
Integer. The number of rows in the x array, as it was declared in the calling program (the leading dimension of x). ldx >= MAX(n1, 1).
ldy
Integer.

The number of rows in the y array, as it was declared in the calling program (the leading dimension of y). ldy >= MAX(n1, 1).

isys
Algorithm used; value dependent on hardware system. Currently, no special options are supported; therefore, you must always specify an isys argument as constant 0.
Output parameters:
y
Array of dimension (0:ldy-1, 0:n2-1). ccfft2d: complex array. zzfft2d: 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 (the input array is overwritten with the transformed values). In this case, it is necessary that ldx = ldy.
table
Real array; dimension (30+ 2 * n1) + (30 + 2 * n2).

Table of factors and trigonometric functions.

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).

work
Real array; dimension 2 * (n1*n2).

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.

static void scfft2d(Int isign, Int n1, Int n2, Float scale, Float* x, Int ldx, Complex* y, Int ldy, Float* table, Float* work, Int isys)
static void dzfft2d(Int isign, Int n1, Int n2, Double scale, Double* x, Int ldx, DComplex* y, Int ldy, Double* table, Double* work, Int isys)
static void csfft2d(Int isign, Int n1, Int n2, Float scale, Complex* x, Int ldx, Float* y, Int ldy, Float* table, Float* work, Int isys)
static void zdfft2d(Int isign, Int n1, Int n2, Double scale, DComplex* x, Int ldx, Double* y, Int ldy, Double* table, Double* work, Int isys)

scfft2d/dzfft2d computes the two-dimensional Fast Fourier Transform (FFT) of the real matrix x, and it stores the results in the complex matrix y. csfft2d/zdfft2d computes the corresponding inverse transform.

In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. First the function of scfft2d is described. Suppose the arrays are dimensioned as follows:

         REAL    X(0:ldx-1, 0:n2-1)
         COMPLEX Y(0:ldy-1, 0:n2-1)
    
    where ldx >= n1 ldy >= (n1/2) + 1.
    

scfft2d computes the formula:

                         n2-1 n1-1
     Y(k1, k2) = scale * Sum  Sum [ X(j1, j2)*w1**(j1*k1)*w2**(j2*k2) ]
                         j2=0 j1=0

     for k1 = 0, ..., n1/2 + 1
         k2 = 0, ..., n2-1

     where:
         w1 = exp(isign*2*pi*i/n1)
         w2 = exp(isign*2*pi*i/n2)
         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 these routines 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 * 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.0/(n1 . n2).

scfft2d is very similar in function to ccfft2d, but it takes the real-to-complex transform in the first dimension, followed by the complex-to-complex transform in the second dimension.

csfft2d does the reverse. It takes the complex-to-complex FFT in the second dimension, followed by the complex-to-real FFT in the first dimension.

See the scfft man page for more information about real-to-complex and complex-to-real FFTs. The two-dimensional analog of the conjugate formula is as follows:

         Y       = conjg Y
          k , k            n1 - k , n2 - k
           1   2                 1        2
    
         for n1/2 <  k  <= n1 - 1
                      1
    
            0 <= k  <= n2 - 1
                  2
    where the notation conjg(z) represents the complex conjugate of z.
    

Thus, you have to compute only (slightly more than) half of the output values, namely:

      Y       for 0 <= k  <= n1/2    0 <= k  <= n2 - 1
       k , k            1                  2
        1   2

Algorithm

scfft2d uses a routine similar to scfftm to do a real-to-complex FFT on the columns, then uses a routine similar to ccfftm to do a complex-to-complex FFT on the rows.

csfft2d uses a routine similar to ccfftm to do a complex-to-complex FFT on the rows, then uses a routine similar to csfftm to do a complex-to-real FFT on the columns.

Table Initialization

The table array stores factors of n1 and n2, and trigonometric tables that are used in calculation of the FFT. table must be initialized by calling the routine with isign = 0. table does not have to be reinitialized if the values of the problem sizes, n1 and n2, do not change.

Dimensions

In the preceding description, it is assumed that array subscripts were zero-based, as is customary in FFT applications. Thus, the input and output arrays are declared:

      REAL    X(0:ldx-1, 0:n2-1)
      COMPLEX Y(0:ldy-1, 0:n2-1)

No change is made in the calling sequence, however, if you prefer to use the more customary Fortran style with subscripts starting at 1. The same values of ldx and ldy would be passed to the subroutine even if the input and output arrays were dimensioned as follows:

      REAL    X(ldx, n2)
      COMPLEX Y(ldy, n2)

Example

The following examples are for Origin series only.

Example 1: Initialize the TABLE array in preparation for doing a two-dimensional FFT of size 128 by 256. In this case, only the isign, n1, n2, and table arguments are used; you can use dummy arguments or zeros for other arguments.

       REAL TABLE ((15 + 128) + 2(15 + 256))
       CALL SCFFT2D (0, 128, 256, 0.0, DUMMY, 1, DUMMY, 1,
      &  TABLE, DUMMY, 0)

Example 2: X is a real array of size (0:128, 0: 255), and Y is a complex array of dimension (0:64, 0:255). The first 128 elements of each column of X contain data; for performance reasons, the extra element forces the leading dimension to be an odd number. Take the two-dimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.

       REAL    X(0:128, 0:255)
       COMPLEX Y(0:64, 0:255)
       REAL    TABLE ((15 + 128) + 2(15 + 256))
       REAL    WORK(128*256)
       ...
       CALL SCFFT2D(0, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
       CALL SCFFT2D(1, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)

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/(128*256) is used. Assume that the TABLE array is initialized already.

       CALL CSFFT2D(-1, 128, 256, 1.0/(128.0*256.0), Y, 65,
      &  X, 130, TABLE, WORK, 0)

Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. No change is needed in the subroutine calls.

       REAL    X(129, 256)
       COMPLEX Y(65, 256)
       ...
       CALL SCFFT2D(0, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)
       CALL SCFFT2D(1, 128, 256, 1.0, X, 129, Y, 65, TABLE, WORK, 0)

Example 5: Perform the same computation as in example 4, but equivalence the input and output arrays to save storage space. In this case, a row must be added to X, because it is equivalenced to a complex array. Assume that TABLE is already initialized.

       REAL    X(130, 256)
       COMPLEX Y(65, 256)
       EQUIVALENCE ( X(1, 1), Y(1, 1) )
       ...
       CALL SCFFT2D(1, 128, 256, 1.0, X, 130, Y, 65, TABLE, WORK, 0)

Input parameters:

isign
Integer. 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, n, and table.

If isign = +1 or -1, the value of isign is the sign of the exponent used in the FFT formula.

n1
Integer. Transform size in the first dimension. If n1 is not positive, scfft2d returns without calculating a transform.
n2
Integer. Transform size in the second dimension. If n2 is not positive, scfft2d returns without calculating a transform.
scale
Scale factor. scfft2d: real. dzfft2d: double precision. csfft2d: real. zdfft2d: double precision. Each element of the output array is multiplied by scale factor after taking the Fourier transform, as defined previously.
x
Array of dimension (0:ldx-1, 0:n2-1). scfft2d: real array. dzfft2d: double precision array. csfft2d: complex array. zdfft2d: double complex array.

Array of values to be transformed.

ldx
Integer. The number of rows in the x array, as it was declared in the calling program. That is, the leading dimension of x. scfft2d, dzfft2d: ldx >= MAX(n1, 1). csfft2d, zdfft2d: ldx >= MAX(n1/2 + 1, 1).
ldy
Integer.

The number of rows in the y array, as it was declared in the calling program (the leading dimension of y).

scfft2d, dzfft2d: ldy >= MAX(n1/2 + 1, 1). csfft2d, zdfft2d: ldy >= MAX(n1 + 2, 1).

In the complex-to-real routine, two extra elements are in the first dimension (ldy >= n1 + 2, rather than just ldy >= n1). These elements are needed for intermediate storage during the computation. On exit, their value is undefined.

isys
Integer array of dimension (0:isys(0)). The first element of the array specifies how many more elements are in the array. Use isys to specify certain processor-specific parameters or options.

If isys(0) = 0, the default values of such parameters are used. In this case, you can specify the argument value as the scalar integer constant 0.

If isys(0) > 0, isys(0) gives the upper bound of the isys array. Therefore, if il = isys(0), user-specified parameters are expected in isys(1) through isys(il).

isys
Algorithm used; value dependent on hardware system. Currently, no special options are supported; therefore, you must always specify an isys argument as constant 0.
Output parameters:
y
scfft2d: complex array. dzfft2d: double complex array. csfft2d: real array. zdfft2d: double precision array.

Output array of transformed values. The output array can be the same as the input array, in which case, the transform is done in place and the input array is overwritten with the transformed values. In this case, it is necessary that the following equalities hold:

scfft2d, dzfft2d: ldx = 2 * ldy. csfft2d, zdfft2d: ldy = 2 * ldx.

table
Real array; dimension (15 + n1) + 2(15 + n2).

Table of factors and trigonometric functions.

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).

work
Real array; dimension (n1 * n2).

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.

static void ccfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Complex* x, Int ldx, Int ldx2, Complex* y, Int ldy, Int ldy2, Float* table, Float* work, Int isys)
static void zzfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, DComplex* x, Int ldx, Int ldx2, DComplex* y, Int ldy, Int ldy2, Double* table, Double* work, Int isys)

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 dimensioned as follows:

      COMPLEX X(0:n1-1, 0:n2-1, 0:n3-1)
      COMPLEX Y(0:n1-1, 0:n2-1, 0:n3-1)

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),
    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).

Example

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.

       REAL TABLE ((30 + 256) + (30 + 256) + (30 + 256))
       CALL CCFFT3D (0, 128, 128, 128, 0.0, DUMMY, 1, 1, DUMMY, 1, 1,
      &  TABLE, DUMMY, 0)

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.

       COMPLEX X(0:128, 0:128, 0:128)
       COMPLEX Y(0:128, 0:128, 0:128)
       REAL TABLE ((30+256) + (30 + 256) + (30 + 256))
       REAL     WORK 2(128*128*128)
       ...
       CALL CCFFT3D(0, 128, 128, 128, 1.0, DUMMY, 1, 1,
     &    DUMMY, 1, 1, TABLE, WORK, 0)
       CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
     &    Y, 129, 129, TABLE, WORK, 0)

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**3) is used. Assume that the TABLE array is already initialized.

       CALL CCFFT3D(-1, 128, 128, 128, 1.0/(128.0**3), Y, 129, 129,
      &  X, 129, 129, TABLE, WORK, 0)

Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. The subroutine calls do not change.

       COMPLEX X(129, 129, 129)
       COMPLEX Y(129, 129, 129)
       ...
       CALL CCFFT3D(0, 128, 128, 128, 1.0, DUMMY, 1, 1,
      &    DUMMY, 1, 1, TABLE, WORK, 0)
       CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
      &    Y, 129, 129, TABLE, WORK, 0)

Example 5: Perform the same computation as in example 4, but put the output back in the array X to save storage space. Assume that the TABLE array is already initialized.

       COMPLEX X(129, 129, 129)
       ...
       CALL CCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
      &    X, 129, 129, TABLE, WORK, 0)

Input parameters:

isign
Integer. 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. Transform size in the first dimension. If n1 is not positive, the routine returns without computing a transform.

n2
Integer. Transform size in the second dimension. If n2 is not positive, the routine returns without computing a transform.

n3
Integer. Transform size in the third dimension. If n3 is not positive, the routine returns without computing a transform.

scale
Scale factor. ccfft3d: real. zzfft3d: double precision.

Each element of the output array is multiplied by scale after taking the Fourier transform, as defined previously.

x
Array of dimension (0:ldx-1, 0:ldx2-1, 0:n3-1). ccfft3d: complex array. zzfft3d: double complex array.

Input array of values to be transformed.

ldx
Integer. The first dimension of x, as it was declared in the calling program (the leading dimension of x). ldx >= MAX(n1, 1).

ldx2
Integer. The second dimension of x, as it was declared in the calling program. ldx2 >= MAX(n2, 1).

ldy
Integer. The first dimension of y, as it was declared in the calling program (the leading dimension of y). ldy >= MAX(n1, 1).

ldy2
Integer. The second dimension of y, as it was declared in the calling program. ldy2 >= MAX(n2, 1).

isys
Algorithm used; value dependent on hardware system. Currently, no special options are supported; therefore, you must always specify an isys argument as constant 0.

isys = 0 or 1 depending on the amount of workspace the user can provide to the routine.

Output parameters:
y
Array of dimension (0:ldy-1, 0:ldy2-1, 0:n3-1). ccfft3d: 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.

table
Real array; dimension 30 + 2 * n1) + (30 + 2 * n2) + (30 + 2 * n3).

Table of factors and trigonometric functions. 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).

work
Real array; dimension (n1 * n2 * n3).

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.

static void scfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Float* x, Int ldx, Int ldx2, Complex* y, Int ldy, Int ldy2, Float* table, Float* work, Int isys)
static void dzfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, Double* x, Int ldx, Int ldx2, DComplex* y, Int ldy, Int ldy2, Double* table, Double* work, Int isys)
static void csfft3d(Int isign, Int n1, Int n2, Int n3, Float scale, Complex* x, Int ldx, Int ldx2, Float* y, Int ldy, Int ldy2, Float* table, Float* work, Int isys)
static void zdfft3d(Int isign, Int n1, Int n2, Int n3, Double scale, DComplex* x, Int ldx, Int ldx2, Double* y, Int ldy, Int ldy2, Double* table, Double* work, Int isys)

These are C++ wrapper functions for the 3D real-to-complex and complex-to-real transform routines in the SGI/Cray Scientific Library (SCSL). The purpose of these definitions is to overload the functions so that C++ users can access the functions in SCSL with identical function names.

Warning Currently, the SCSL is available only on SGI machines.

scfft3d/dzfft3d computes the three-dimensional Fast Fourier Transform (FFT) of the real matrix X, and it stores the results in the complex matrix Y. csfft3d/zdfft3d computes the corresponding inverse transform.

In FFT applications, it is customary to use zero-based subscripts; the formulas are simpler that way. First, the function of SCFFT3D is described. Suppose the arrays are dimensioned as follows:

      REAL    X(0:ldx-1, 0:ldx2-1, 0:n3-1)
      COMPLEX Y(0:ldy-1, 0:ldy2-1, 0:n3-1)

scfft3d computes 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/2,
        k2 = 0, ..., n2 - 1,
        k3 = 0, ..., n3 - 1,

    where:
        w1 = exp(isign*2*pi*i/n1),
        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 these routines 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 isign = -1 and

      scale = 1.0/(n1*n2*n3).

scfft3d is very similar in function to ccfft3d, but it takes the real-to-complex transform in the first dimension, followed by the complex-to-complex transform in the second and third dimensions.

csfft3d does the reverse. It takes the complex-to-complex FFT in the third and second dimensions, followed by the complex-to-real FFT in the first dimension.

See the scfftm man page for more information about real-to-complex and complex-to-real FFTs. The three dimensional analog of the conjugate formula is as follows:

         Y         = conjg Y
          k ,k ,k            n1 - k , n2 - k , n3 - k
           1  2  3                 1        2        3
    
         for  n1/2 <  k  <= n1 - 1
                       1
    
              0 <= k  <= n2 - 1
                    2
    
              0 <= k  <= n3 - 1
                    3
    where the notation conjg(z) represents the complex conjugate of z.
    

Thus, you have to compute only (slightly more than) half out the output values, namely:

      Y
       k ,k ,k
        1  2  3

      for  0 <= k  <= n1/2
                 1

           0 <= k  <= n2 - 1
                 2

           0 <= k  <= n3 - 1

Algorithm

scfft3d uses a routine similar to scfftm to do multiple FFTs first on all columns of the input matrix, then uses a routine similar to ccfftm on all rows of the result, and then on all planes of that result. See scfftm and ccfftm for more information about the algorithms used.

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.

       REAL TABLE ((15 + 128) + 2(15+128) + 2( 15 + 128))
       CALL SCFFT3D (0, 128, 128, 128, 0.0, DUMMY, 1, 1, DUMMY, 1, 1,
      &  TABLE, DUMMY, 0)

Example 2: X is a real array of size (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. Y is a complex array of dimension (0:64, 0:128, 0:128). Take the three-dimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.

       REAL    X(0:128, 0:128, 0:128)
       COMPLEX Y(0:64,  0:128, 0:128)
       REAL TABLE ((15+128) + 2(15 + 128) + 2(15 + 128))
       REAL    WORK(128*128*128)
       ...
       CALL SCFFT3D(0, 128, 128, 128, 1.0, X, 129, 129,
      &    Y, 65, 129, TABLE, WORK, 0)
       CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
      &    Y, 65, 129, TABLE, WORK, 0)

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/(128**3) is used. Assume that the TABLE array is initialized already.

       CALL CSFFT3D(-1, 128, 128, 128, 1.0/128.0**3, Y, 65, 129,
      &    X, 130, 129, TABLE, WORK, 0)

Example 4: Perform the same computation as in example 2, but assume that the lower bound of each array is 1, rather than 0. No change is made in the subroutine calls.

       REAL    X(129, 129, 129)
       COMPLEX Y(65, 129, 129)
       REAL TABLE ((15+128) + 2(15 + 128) + 2(15 + 128))
       REAL    WORK(128*128*128)
       ...
       CALL SCFFT3D(0, 128, 128, 128, 1.0, X, 129, 129,
      &    Y, 65, 129, TABLE, WORK, 0)
       CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 129, 129,
      &    X, 129, 129, TABLE, WORK, 0)

Example 5: Perform the same computation as in example 4, but equivalence the input and output arrays to save storage space. Assume that the TABLE array is initialized already.

       REAL    X(130, 129, 129)
       COMPLEX Y(65, 129, 129)
       EQUIVALENCE (X(1, 1, 1), Y(1, 1, 1))
       ...
       CALL SCFFT3D(1, 128, 128, 128, 1.0, X, 130, 129,
      &    Y, 65, 129, TABLE, WORK, 0)

Input parameters:

isign
Integer. 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. Transform size in the first dimension. If n1 is not positive, scfft3d returns without computing a transform.

n2
Integer. Transform size in the second dimension. If n2 is not positive, scfft3d returns without computing a transform.

n3
Integer. Transform size in the third dimension. If n3 is not positive, scfft3d returns without computing a transform.

scale
Scale factor. scfft3d: real. dzfft3d: double precision. csfft3d: real. zdfft3d: double precision. Each element of the output array is multiplied by scale after taking the Fourier transform, as defined previously.

x
Array of dimension (0:ldx-1, 0:ldx2-1, 0:n3-1). scfft3d: real array. dzfft3d: double precision array. csfft3d: complex array. zdfft3d: double complex array.

Array of values to be transformed.

ldx
Integer. The first dimension of x, as it was declared in the calling program (the leading dimension of x).

scfft3d, dzfft3d: ldx >= MAX(n1, 1). csfft3d, zdfft3d: ldx >= MAX(n1/2 + 1, 1).

ldx2
Integer. The second dimension of x, as it was declared in the calling program. ldx2 >= MAX(n2, 1).

ldy
Integer. The first dimension of y, as it was declared in the calling program; that is, the leading dimension of y.

scfft3d, dzfft3d: ldy >= MAX(n1/2 + 1, 1). csfft3d, zdfft3d: ldy >= MAX(n1 + 2, 1).

In the complex-to-real routine, two extra elements are in the first dimension (that is, ldy >= n1 + 2, rather than just ldy >= n1). These elements are needed for intermediate storage during the computation. On exit, their value is undefined.

ldy2
Integer. The second dimension of y, as it was declared in the calling program. ldy2 >= MAX(n2, 1).

isys
Algorithm used; value dependent on hardware system. Currently, no special options are supported; therefore, you must always specify an isys argument as constant 0.

isys = 0 or 1 depending on the amount of workspace the user can provide to the routine.

Output parameters:
y
Array of dimension (0:ldy-1, 0:ldy2-1, 0:n3-1). scfft3d: complex array. dzfft3d: double complex array. csfft3d: real array. zdfft3d: double precision array.

Output array of transformed values. The output array can 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 the following equalities hold:

scfft3d, dzfft3d: ldx = 2 * ldy, and ldx2 = ldy2. csfft3d, zdfft3d: ldy = 2 * ldx, and ldx2 = ldy2.

table
Real array; dimension (15 + n1) + 2(15 + n2) + 2(15 + n3).

Table of factors and trigonometric functions.

This array must be initialized by a call to scfft3d or csfft3d with isign = 0.

If isign = 0, table is initialized to contain trigonometric tables needed to compute a three-dimensional FFT of size n1 by n2 by n3. If isign = +1 or -1, the values in table are assumed to be initialized already by a prior call with isign = 0.

work
Real array; dimension n1 * n2 * n3.

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.