SCRFT and DCRFT (Complex-to-Real Fourier Transform)
Purpose
These subroutines compute a set of m real discrete n-point Fourier transforms of complex conjugate even data.
Data Types | ||
---|---|---|
X |
Y , scale |
Subroutine |
Short-precision complex | Short-precision real | SCRFT |
Long-precision complex | Long-precision real | DCRFT |
- Two invocations of this subroutine are necessary: one to prepare the working storage for the subroutine, and the other to perform the computations.
- On certain processors, SIMD algorithms may be used if alignment requirements are met. For further details, see Use of SIMD Algorithms by Some Subroutines in the Libraries Provided by ESSL.
Syntax
Language | Syntax |
---|---|
Fortran | CALL SCRFT (init, x, inc2x,
y, inc2y, n, m,
isign, scale, aux1,
naux1, aux2, naux2,
aux3, naux3) CALL DCRFT (init, x, inc2x, y, inc2y, n, m, isign, scale, aux1, naux1, aux2, naux2) |
C and C++ | scrft (init, x, inc2x,
y, inc2y, n, m,
isign, scale, aux1,
naux1, aux2, naux2,
aux3, naux3); dcrft (init, x, inc2x, y, inc2y, n, m, isign, scale, aux1, naux1, aux2, naux2); |
- On Entry
- init
- is a flag, where:
If init ≠ 0, trigonometric functions and other parameters, depending on arguments other than x, are computed and saved in aux1. The contents of x and y are not used or changed.
If init = 0, the discrete Fourier transforms of the given sequences are computed. The only arguments that may change after initialization are x, y, and aux2. All scalar arguments must be the same as when the subroutine was called for initialization with init ≠ 0.
Specified as: an integer. It can have any value.
- x
- is the array
X
, consisting of m sequences. Due to complex conjugate symmetry, the input consists of only the first (n/2)+1 elements of each sequence; that is, xji, j = 0, 1, …, n/2, i = 1, 2, …, m. The sequences are assumed to be stored with stride 1.Specified as: an array of (at least) length n/2+1+(m-1)inc2x, containing numbers of the data type indicated in Table 1. This array can be declared asX
(inc2x,m). - inc2x
- is the stride between the first elements of the sequences in array
X
. (If m = 1, this argument is ignored.) Specified as: an integer; inc2x ≥ (n/2)+1. - y
- See On Return.
- inc2y
- is the stride between the first elements of the sequences in array
Y
. (If m = 1, this argument is ignored.) Specified as: an integer; inc2y ≥ n. - n
- is the length of each sequence to be transformed.
Specified as: an integer; n ≤ 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.
- m
- is the number of sequences to be transformed.
Specified as: an integer; m > 0.
- isign
- controls the direction of the transform, determining the sign Isign of the
exponent of Wn, where:
If isign = positive value, Isign = + (transforming time to frequency).
If isign = negative value, Isign = - (transforming frequency to time).
Specified as: an integer; isign > 0 or isign < 0.
- scale
- is the scaling constant scale. See Function for its usage.
- aux1
- is the working storage for this subroutine, where:
If init ≠ 0, the working storage is computed.
If init = 0, the working storage is used in the computation of the Fourier transforms.
Specified as: an area of storage, containing naux1 long-precision real numbers.
- naux1
- is the number of doublewords in the working storage specified in aux1.
Specified as: an integer; naux1 > 13 (32-bit integer arguments) or 25 (64-bit integer arguments) and naux1 ≥ (minimum value required for successful processing). To determine a sufficient value, use the processor-independent formulas. For values between 13 (32-bit integer arguments) or 25 (64-bit integer arguments) and the minimum value, you have the option of having the minimum value returned in this argument. For details, see Using Auxiliary Storage in ESSL.
- aux2
- has the following meaning:
If naux2 = 0 and error 2015 is unrecoverable, aux2 is ignored.
Otherwise, it is the working storage used by this subroutine that is available for use by the calling program between calls to this subroutine.
Specified as: an area of storage, containing naux2 long-precision real numbers. On output, the contents are overwritten.
- naux2
- is the number of doublewords in the working storage specified in aux2.
Specified as: an integer, where:
If naux2 = 0 and error 2015 is unrecoverable, SCRFT and DCRFT dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.
Otherwise, naux2 ≥ (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.
- aux3
- this argument is provided for migration purposes only and is ignored.
Specified as: an area of storage, containing naux3 long-precision real numbers.
- naux3
- this argument is provided for migration purposes only and is ignored.
Specified as: an integer.
- On Return
- y
- has the following meaning, where:
If init ≠ 0, this argument is not used, and its contents remain unchanged.
If init = 0, this is array
Y
, consisting of the results of the m discrete Fourier transforms of the complex conjugate even data, each of length n. The sequences are stored with stride 1. - aux1
- is the working storage for this subroutine, where:
If init ≠ 0, it contains information ready to be passed in a subsequent invocation of this subroutine.
If init = 0, its contents are unchanged.
Returned as: the contents are not relevant.
Notes
- aux1 should not be used by the calling program between calls to this subroutine with init ≠ 0 and init = 0. However, it can be reused after intervening calls to this subroutine with different arguments.
- When using the ESSL SMP Libraries, for optimal performance, the number of threads specified should be the same for init ≠ 0 and init = 0.
- The elements in each sequence in x and y are assumed to be stored in
contiguous storage locations—that
is, with a stride of 1. Therefore, inc1x and inc1y values are
not a part of the argument list. For optimal performance, the inc2x and
inc2y values should be close to their respective minimum values, which are given
below: min(inc2y) = n
min(inc2x) = n/2+1If you specify the same array for
X
andY
, then inc2y must equal 2(inc2x). In this case, output overwrites input. If m = 1, the inc2x and inc2y values are not used by the subroutine. If you specify different arrays forX
andY
, they must have no common elements; otherwise, results are unpredictable. See Vector concepts.
Formulas
Processor-Independent Formulas for SCRFT for NAUX1 and NAUX2:
- NAUX1 Formulas
-
- For 32-bit integer arguments:
-
If n ≤ 16384, use naux1 = 25000.
If n > 16384, use naux1 = 20000+0.82n. - For 64-bit integer arguments:
-
If n ≤ 16384, use naux1 = 35000.
If n > 16384, use naux1 = 30000+0.82n.
- NAUX2 Formulas
-
If n ≤ 16384, use naux2 = 20000.
If n > 16384, use naux2 = 20000+0.57n.
Processor-Independent Formulas for DCRFT for NAUX1 and NAUX2:
- NAUX1 Formulas
-
- For 32-bit integer arguments:
-
If n ≤ 4096, use naux1 = 22000.
If n > 4096, use naux1 = 20000+1.64n. - For 64-bit integer arguments:
-
If n ≤ 4096, use naux1 = 32000.
If n > 4096, use naux1 = 30000+1.64n.
- NAUX2 Formulas
-
If n ≤ 4096, use naux2 = 20000.
If n > 4096, use naux2 = 20000+1.14n.
Function
The set of m real discrete n-point Fourier transforms of
complex conjugate even data in array X
, with results going into array
Y
, is expressed as follows:
for:
i = 1, 2, …, m
where:
and where:
X
.yki are elements of the sequences in array
Y
.Isign is + or - (determined by argument isign).
scale is a scalar value.
Because of the symmetry, Y
has real data. 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/n and isign being
negative. See references [1], [5], [27], and [28].
- With init ≠ 0, the subroutine tests and initializes arguments of the program, setting up the aux1 working storage.
- With init = 0, the subroutine checks that the initialization arguments in the aux1 working storage correspond to the present arguments, and if so, performs the calculation of the Fourier transforms.
Error conditions
- Resource Errors
- Error 2015 is unrecoverable, naux2 = 0, and unable to allocate work area.
- Computational Errors
- None
- Input-Argument Errors
-
- n > 37748736
- m ≤ 0
- inc2x < n/2+1
- inc2y < n
- scale = 0.0
- isign = 0
- The subroutine has not been initialized with the present arguments.
- The length of the transform in n is not an allowable value. Return code 1 is returned if error 2030 is recoverable.
- naux1 ≤ 13
- naux1 is too small—that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.
- Error 2015 is recoverable or naux2≠0, and naux2 is too small—that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.
Examples
- Example 1
-
This example uses the mixed-radix capability and shows how to compute a single transform. The arrays are declared as follows:
COMPLEX*8 X(0:6) REAL*8 AUX1(50), AUX2(1), AUX3(1) REAL*4 Y(0:11)
First, initializeAUX1
using the calling sequence shown below withINIT
≠ 0. Then use the same calling sequence withINIT
= 0 to do the calculation.Note:X
shows the n/2+1 = 7 elements used in the computation.- Because
NAUX2
= 0, this subroutine dynamically allocates theAUX2
working storage.
Call Statement and Input:INIT X INC2X Y INC2Y N M ISIGN SCALE AUX1 NAUX1 AUX2 NAUX2 AUX3 NAUX3 | | | | | | | | | | | | | | | CALL SCRFT(INIT, X , 0 , Y , 0 , 12 , 1 , 1 , SCALE, AUX1 , 50 , AUX2 , 0 , AUX3 , 0 )
INIT = 1
(for initialization)
INIT = 0
(for computation)
SCALE = 1.0
X
contains the following sequence:(1.0, 0.0) (0.0, 0.0) (0.0, 0.0) (0.0, 0.0) (0.0, 0.0) (0.0, 0.0) (0.0, 0.0)
Output:Y = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
- Example 2
-
This example shows another transform computation with different data using the same initialized array
AUX1
as in Example 1.Note: BecauseNAUX2
= 0, this subroutine dynamically allocates theAUX2
working storage.Call Statement and Input:INIT X INC2X Y INC2Y N M ISIGN SCALE AUX1 NAUX1 AUX2 NAUX2 AUX3 NAUX3 | | | | | | | | | | | | | | | CALL SCRFT( 0 , X , 0 , Y , 0 , 12 , 1 , 1 , SCALE, AUX1 , 50 , AUX2 , 0 , AUX3 , 0 )
SCALE = 1.0
X
contains the following sequence:(1.0, 0.0) (1.0, 0.0) (1.0, 0.0) (1.0, 0.0) (1.0, 0.0) (1.0, 0.0) (1.0, 0.0)
Output:Y = (12.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0)
- Example 3
-
This example shows how to compute many transforms simultaneously. The arrays are declared as follows:
COMPLEX*8 X(0:8,2) REAL*8 AUX1(50), AUX2(1), AUX3(1) REAL*4 Y(0:15,2)
First, initialize
AUX1
using the calling sequence shown below withINIT
≠ 0. Then use the same calling sequence withINIT
= 0 to do the calculation.Note: BecauseNAUX2
= 0, this subroutine dynamically allocates theAUX2
working storage.Call Statement and Input:INIT X INC2X Y INC2Y N M ISIGN SCALE AUX1 NAUX1 AUX2 NAUX2 AUX3 NAUX3 | | | | | | | | | | | | | | | CALL SCRFT(INIT, X , 9 , Y , 16 , 16 , 2 , 1 , SCALE, AUX1 , 50 , AUX2 , 0 , AUX3 , 0 )
INIT = 1
(for initialization)
INIT = 0
(for computation)
SCALE = 1.0
X
contains the following two sequences:(1.0, 0.0) (0.0, 0.0) (1.0, 0.0) (0.0, 0.0) (1.0, 0.0) (0.0, 0.0) (1.0, 0.0) (0.0, 0.0) (1.0, 0.0) (0.0, 0.0) (1.0, 0.0) (0.0, 0.0) (1.0, 0.0) (0.0, 0.0) (1.0, 0.0) (0.0, 0.0) (1.0, 0.0) (1.0, 0.0)
Output:
Y
contains the following two sequences:16.0 1.0 0.0 -1.0 0.0 1.0 0.0 -1.0 0.0 1.0 0.0 -1.0 0.0 1.0 0.0 -1.0 0.0 1.0 0.0 -1.0 0.0 1.0 0.0 -1.0 0.0 1.0 0.0 -1.0 0.0 1.0 0.0 -1.0
- Example 4
-
This example shows the same array being used for input and output. The arrays are declared as follows:
COMPLEX*16 X(0:8,2) REAL*8 AUX1(50), AUX2(1) REAL*8 Y(0:17,2)
ArraysX
andY
are made equivalent by the following statement, making them occupy the same storage:EQUIVALENCE (X,Y)
This requires
INC2Y
= 2(INC2X
). First, initializeAUX1
using the calling sequence shown below withINIT
≠ 0. Then use the same calling sequence withINIT
= 0 to do the calculation.Note: BecauseNAUX2
= 0, this subroutine dynamically allocates theAUX2
working storage.Call Statement and Input:INIT X INC2X Y INC2Y N M ISIGN SCALE AUX1 NAUX1 AUX2 NAUX2 | | | | | | | | | | | | | CALL DCRFT(INIT, X , 9 , Y , 18 , 16 , 2 , -1 , SCALE, AUX1 , 50 , AUX2 , 0)
INIT = 1
(for initialization)
INIT = 0
(for computation)
SCALE = 0.0625
X
contains the following two sequences:( 1.0, 0.0) ( 1.0, 0.0) ( 0.0, 1.0) ( 0.0, -1.0) (-1.0, 0.0) (-1.0, 0.0) ( 0.0, -1.0) ( 0.0, 1.0) ( 1.0, 0.0) ( 1.0, 0.0) ( 0.0, 1.0) ( 0.0, -1.0) (-1.0, 0.0) (-1.0, 0.0) ( 0.0, -1.0) ( 0.0, 1.0) ( 1.0, 0.0) ( 1.0, 0.0)
Output:
Y
contains the following two sequences:0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0