Engineering and Scientific Subroutine Library for AIX Version 3 Release 3: Guide and Reference

SCFT3 and DCFT3--Complex Fourier Transform in Three Dimensions

These subroutines compute the three-dimensional discrete Fourier transform of complex data.

Table 139. Data Types

Note:

For each use, only one invocation of this subroutine is necessary. The initialization phase, preparing the working storage, is a relatively small part of the total computation, so it is performed on each invocation.

Syntax

On Entry

x

is the array X, containing the three-dimensional data to be transformed, where each element x_j1,j2,j3, using zero-based indexing, is stored in X(j1+j2(inc2x)+j3(inc3x)) for j1 = 0, 1, ..., n1-1, j2 = 0, 1, ..., n2-1, and j3 = 0, 1, ..., n3-1. The strides for the elements in the first, second, and third dimensions are assumed to be 1, inc2x( >= n1), and inc3x( >= (n2)(inc2x)), respectively.

Specified as: an array, containing numbers of the data type indicated in Table 139. This array must be aligned on a doubleword boundary. If the array is dimensioned X(LDA1,LDA2,LDA3), then LDA1 = inc2x, (LDA1)(LDA2) = inc3x, and LDA3 >= n3. For information on how to set up this array, see Setting Up Your Data. For more details, see Notes.

inc2x

is the stride between the elements in array X for the second dimension. Specified as: a fullword integer; inc2x >= n1.

inc3x

is the stride between the elements in array X for the third dimension. Specified as: a fullword integer; inc3x >= (n2)(inc2x).

y

See On Return.

inc2y

is the stride between the elements in array Y for the second dimension. Specified as: a fullword integer; inc2y >= n1.

inc3y

is the stride between the elements in array Y for the third dimension. Specified as: a fullword integer; inc3y >= (n2)(inc2y).

n1

is the length of the first dimension of the three-dimensional data in the array to be transformed. Specified as: a fullword integer; n1 <= 37748736 and must be one of the values listed in Acceptable Lengths for the Transforms. For all other values specified less than 37748736, you have the option of having the next larger acceptable value returned in this argument. For details, see Providing a Correct Transform Length to ESSL.

n2

is the length of the second dimension of the three-dimensional data in the array to be transformed. Specified as: a fullword integer; n2 <= 37748736 and must be one of the values listed in Acceptable Lengths for the Transforms. For all other values specified less than 37748736, you have the option of having the next larger acceptable value returned in this argument. For details, see Providing a Correct Transform Length to ESSL.

n3

is the length of the third dimension of the three-dimensional data in the array to be transformed. Specified as: a fullword integer; n3 <= 37748736 and must be one of the values listed in Acceptable Lengths for the Transforms. For all other values specified less than 37748736, you have the option of having the next larger acceptable value returned in this argument. For details, see Providing a Correct Transform Length to ESSL.

isign

controls the direction of the transform, determining the sign Isign of the exponents of W_n1, W_n2, and W_n3, where:

If isign = positive value, Isign = + (transforming time to frequency).

If isign = negative value, Isign = - (transforming frequency to time).

Specified as: a fullword integer; isign > 0 or isign < 0.

scale

is the scaling constant scale. See Function for its usage. Specified as: a number of the data type indicated in Table 139, where scale > 0.0 or scale < 0.0.

aux

has the following meaning:

If naux = 0 and error 2015 is unrecoverable, aux is ignored.

Otherwise, it is a storage work area used by this subroutine.

Specified as: an area of storage, containing naux long-precision real numbers. On output, the contents are overwritten.

naux

is the number of doublewords in the working storage specified in aux. Specified as: a fullword integer, where:

If naux = 0 and error 2015 is unrecoverable, SCFT3 and DCFT3 dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, naux >=  (minimum value required for successful processing). To determine a sufficient value, use the processor-independent formulas. For all other values specified less than the minimum value, you have the option of having the minimum value returned in this argument. For details, see Using Auxiliary Storage in ESSL.

On Return

y

is the array Y, containing the elements resulting from the three-dimensional discrete Fourier transform of the data in X. Each element y_k1,k2,k3, using zero-based indexing, is stored in Y(k1+k2(inc2y)+k3(inc3y)) for k1 = 0, 1, ..., n1-1, k2 = 0, 1, ..., n2-1, and k3 = 0, 1, ..., n3-1. The strides for the elements in the first, second, and third dimensions are assumed to be 1, inc2y( >= n1), and inc3y( >= (n2)(inc2y)), respectively.

Returned as: an array, containing numbers of the data type indicated in Table 139. This array must be aligned on a doubleword boundary. If the array is dimensioned Y(LDA1,LDA2,LDA3), then LDA1 = inc2y, (LDA1)(LDA2) = inc3y, and LDA3 >= n3. For information on how to set up this array, see Setting Up Your Data. For more details, see Notes.

Notes

1. If you specify the same array for X and Y, then inc2x must be greater than or equal to inc2y, and inc3x must be greater than or equal to inc3y. In this case, output overwrites input. When using the ESSL SMP library in a multithreaded environment, if inc2x > inc2y or inc3x > inc3y, these subroutines run on a single thread and issue an attention message.

   If you specify different arrays X and Y, they must have no common elements; otherwise, results are unpredictable. See Concepts.

2. You should use STRIDE--Determine the Stride Value for Optimal Performance in Specified Fourier Transform Subroutines to determine the optimal values for the strides inc2y and inc3y for your output array. The strides for your input array do not affect performance. Example 7 in the STRIDE subroutine description explains how it is used for these subroutines. For additional information on how to set up your data, see Setting Up Your Data.

Processor-Independent Formulas for SCFT3 for NAUX

Use the following formulas for calculating naux:

1. If max(n2, n3) < 252 and:

   If n1 <= 8192, use naux = 60000.

   If n1 > 8192, use naux = 60000+2.28n1.

2. If n2 >= 252, n3 < 252, and:

   If n1 <= 8192, use naux = 60000+lambda.

   If n1 > 8192, use naux = 60000+2.28n1+lambda,

   where lambda = (n2+256)(s+2.28)

   and s = min(64, n1).

3. If n2 < 252, n3 >= 252, and:

   If n1 <= 8192, use naux = 60000+psi.

   If n1 > 8192, use naux = 60000+2.28n1+psi,

   where psi = (n3+256)(s+2.28)

   and s = min(64, (n1)(n2)).

4. If n2 >= 252 and n3 >= 252, use the larger of the values calculated for cases 2 and 3 above.

Processor-Independent Formulas for DCFT3 for NAUX

Use the following formulas for calculating naux:

1. If max(n2, n3) < 252 and:

   If n1 <= 2048, use naux = 60000.

   If n1 > 2048, use naux = 60000+4.56n1.

2. If n2 >= 252, n3 < 252, and:

   If n1 <= 2048, use naux = 60000+lambda.

   If n1 > 2048, use naux = 60000+4.56n1+lambda,

   where lambda = ((2)n2+256)(s+4.56)

   and s = min(64, n1).

3. If n2 < 252, n3 >= 252, and:

   If n1 <= 2048, use naux = 60000+psi.

   If n1 > 2048, use naux = 60000+4.56n1+psi,

   where psi = ((2)n3+256)(s+4.56)

   and s = min(64, (n1)(n2)).

4. If n2 >= 252 and n3 >= 252, use the larger of the values calculated for cases 2 and 3 above.

Function

The three-dimensional discrete Fourier transform of complex data in array X, with results going into array Y, is expressed as follows:

for:

k1 = 0, 1, ..., n1-1

k2 = 0, 1, ..., n2-1

k3 = 0, 1, ..., n3-1

where:

and where:

x_j1,j2,j3 are elements of array X.

y_k1,k2,k3 are elements of array Y.

Isign is + or - (determined by argument isign).

scale is a scalar value.

For scale = 1.0 and isign being positive, you obtain the discrete Fourier transform, a function of frequency. The inverse Fourier transform is obtained with scale = 1.0/((n1)(n2)(n3)) and isign being negative. See references [1], [4], [5], [19], and [20].

Error Conditions

Resource Errors

Error 2015 is unrecoverable, naux = 0, and unable to allocate work area.

Computational Errors

None

Input-Argument Errors

1. n1 > 37748736
2. n2 > 37748736
3. n3 > 37748736
4. inc2x < n1
5. inc3x < (n2)(inc2x)
6. inc2y < n1
7. inc3y < (n2)(inc2y)
8. scale = 0.0
9. isign = 0
10. The length of one of the transforms in n1, n2, or n3 is not an allowable value. Return code 1 is returned if error 2030 is recoverable.
11. Error 2015 is recoverable or naux<>0, and naux is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.

Example

This example shows how to compute a three-dimensional transform. In this example, INC2X >= INC2Y and INC3X >= INC3Y, so that the same array can be used for both input and output. The STRIDE subroutine is called to select good values for the INC2Y and INC3Y strides. (As explained below, STRIDE is not called for INC2X and INC3X.) Using the transform lengths (N1 = 32, N2 = 64, and N3 = 40) along with the output data type (short-precision complex: 'C'), STRIDE is called once for each stride needed. First, it is called for INC2Y:

   CALL STRIDE (N2,N1,INC2Y,'C',0)

The output value returned for INC2Y is 32. Then STRIDE is called again for INC3Y:

   CALL STRIDE (N3,N2*INC2Y,INC3Y,'C',0)

The output value returned for INC3Y is 2056. Because INC3Y is not a multiple of INC2Y, Y is not declared as a three-dimensional array. It is declared as a two-dimensional array, Y(INC3Y,N3).

To equivalence the X and Y arrays requires INC2X >= INC2Y and INC3X >= INC3Y. Therefore, INC2X is set equal to INC2Y( = 32). Also, to declare the X array as a three-dimensional array, INC3X must be a multiple of INC2X. Therefore, its value is set as INC3X = (65)(INC2X) = 2080.

The arrays are declared as follows:

     COMPLEX*8  X(32,65,40),Y(2056,40)
     REAL*8     AUX(30000)

Arrays X and Y are made equivalent by the following statement, making them occupy the same storage:

     EQUIVALENCE (X,Y)

Call Statement and Input

            X  INC2X INC3X  Y  INC2Y INC3Y  N1   N2   N3 ISIGN SCALE   AUX   NAUX
            |    |     |    |    |     |     |    |    |   |     |      |      |
CALL SCFT3( X , 32 , 2080 , Y , 32 , 2056 , 32 , 64 , 40 , 1 , SCALE , AUX , 30000)

Output

Y has (1.0,2.0) in locations Y(ij,k), where ij = 1, 2048 and j = 1, 40. It remains unchanged elsewhere.