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

SRCFT3 and DRCFT3--Real-to-Complex Fourier Transform in Three Dimensions

These subroutines compute the three-dimensional discrete Fourier transform of real data in a three-dimensional array.

Table 140. 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 140. 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/2+1.

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 140, 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, SRCFT3 and DRCFT3 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/2, k2 = 0, 1, ..., n2-1, and k3 = 0, 1, ..., n3-1. Due to complex conjugate symmetry, the output consists of only the first n1/2+1 values along the first dimension of the array, for k1 = 0, 1, ..., n1/2. The strides for the elements in the first, second, and third dimensions are assumed to be 1, inc2y( >= n1/2+1), and inc3y( >= (n2)(inc2y)), respectively.

Returned as: an array, containing numbers of the data type indicated in Table 140. 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 (2)(inc2y), and inc3x must be greater than or equal to (2)(inc3y). In this case, output overwrites input. When using the ESSL SMP library in a multithreaded environment, if inc2x > (2)(inc2y) or inc3x > (2)(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. To achieve the best performance, you should align array X on a doubleword boundary, and inc2x and inc3x should be even numbers. The strides for your input array do not affect performance as long as they are even numbers. In addition, 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. Example 8 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 SRCFT3 for NAUX

Use the following formulas for calculating naux:

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

   If n1 <= 16384, use naux = 65000.

   If n1 > 16384, use naux = 60000+1.39n1.

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

   If n1 <= 16384, use naux = 65000+lambda.

   If n1 > 16384, use naux = 60000+1.39n1+lambda,

   where lambda = (n2+256)(s+2.28) and s = min(64, 1+n1/2).

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

   If n1 <= 16384, use naux = 65000+psi.

   If n1 > 16384, use naux = 60000+1.39n1+psi,

   where psi = (n3+256)(s+2.28) and s = min(64, (n2)(1+n1/2)).

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

If inc2x or inc3x is an odd number, or if array X is not aligned on a doubleword boundary, you should add the following amount to all the formulas given above:

n2(1+n1/2)

Processor-Independent Formulas for DRCFT3 for NAUX

Use the following formulas for calculating naux:

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

   If n1 <= 4096, use naux = 62000.

   If n1 > 4096, use naux = 60000+2.78n1.

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

   If n1 <= 4096, use naux = 62000+lambda.

   If n1 > 4096, use naux = 60000+2.78n1+lambda,

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

   and s = min(64, n1/2).

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

   If n1 <= 4096, use naux = 62000+psi.

   If n1 > 4096, use naux = 60000+2.78n1+psi,

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

   and s = min(64, n2(1+n1/2)).

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

Function

The three-dimensional complex conjugate even discrete Fourier transform of real 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.

The output in array Y is complex. 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/2+1
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 >= (2)(INC2Y) and INC3X >= (2)(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. Using the transform lengths (N1 = 32, N2 = 64, and N2 = 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/2+1,INC2Y,'C',0)

The output value returned for INC2Y is 17. (This value is equal to N1/2+1.) Then STRIDE is called again for INC3Y:

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

The output value returned for INC3Y is 1088. Because INC3Y is a multiple of INC2Y--that is, INC3Y = (N2)(INC2Y)--Y is declared as a three-dimensional array, Y(17,64,40). (In general, for larger arrays, these types of values for INC2Y and INC3Y are not returned by STRIDE, and you are probably not able to declare Y as a three-dimensional array.)

To equivalence the X and Y arrays requires INC2X >= (2)(INC2Y) and INC3X >= (2)(INC3Y). Therefore, the values INC2X = (2)(INC2Y) = 34 and INC3X = (2)(INC3Y) = 2176 are set, and X is declared as a three-dimensional array, X(34,64,40).

The arrays are declared as follows:

     REAL*4     X(34,64,40)
     COMPLEX*8  Y(17,64,40)
     REAL*8     AUX(32000)

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 SRCFT3( X , 34 , 2176 , Y , 17 , 1088 , 32 , 64 , 40 , 1 , SCALE , AUX , 32000)

Output

Y has (1.0,0.0) in all locations.