casa
5.7.016

C++ Interface to the Sgi/Cray Scientific Library (SCSL) More...
#include <SCSL.h>
Static Public Member Functions  
static void  ccfft (Int isign, Int n, Float scale, Complex *x, Complex *y, Float *table, Float *work, Int isys) 
These routines compute the Fast Fourier Transform (FFT) of the complex vector x, and store the result in vector y. More...  
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) 
scfft/dzfft computes the FFT of the real array x, and it stores the results in the complex array y. More...  
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) 
<tt>ccfftm/zzfftm</tt> computes the FFT of each column of the complex matrix x, and stores the results in the columns of complex matrix y. More...  
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) 
<tt>scfftm/dzfftm</tt> 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. More...  
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) 
These routines compute the twodimensional complex Fast Fourier Transform (FFT) of the complex matrix x, and store the results in the complex matrix y. More...  
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) 
<tt>scfft2d/dzfft2d</tt> computes the twodimensional Fast Fourier Transform (FFT) of the real matrix x, and it stores the results in the complex matrix y. More...  
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) 
These routines compute the threedimensional complex FFT of the complex matrix X, and store the results in the complex matrix Y. More...  
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) 
These are C++ wrapper functions for the 3D realtocomplex and complextoreal transform routines in the SGI/Cray Scientific Library (SCSL). More...  
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) 
C++ Interface to the Sgi/Cray Scientific Library (SCSL)
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;

static 
These routines compute the Fast Fourier Transform (FFT) of the complex vector x, and store the result in vector y.
ccfft
does the complextocomplex transform and zzfft
does the same for double precision arrays.
In FFT applications, it is customary to use zerobased subscripts; the formulas are simpler that way. Suppose that the arrays are dimensioned as follows:
The output array is the FFT of the input array, using the following formula for the FFT:
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.
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.
In the preceding description, it is assumed that array subscripts were zerobased, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:
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):
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.
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.
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.
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.
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.
Input parameters:
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.
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. Array of dimension (0:n1). ccfft
: complex array. zzfft
: double complex array.
Input array of values to be transformed.
Output parameters:
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. 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).
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 

static 
These routines compute the twodimensional complex Fast Fourier
Transform (FFT) of the complex matrix x, and store the results in the complex matrix y.
ccfft2d
does the complextocomplex transform and zzfft
does the same for double precision arrays.
In FFT applications, it is customary to use zerobased subscripts; the formulas are simpler that way. Suppose that the arrays are dimensioned as follows:
These routines compute the formula:
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).
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.
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.
In the preceding description, it is assumed that array subscripts were zerobased, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:
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:
All examples here are for Origin series only.
Example 1: Initialize the TABLE array in preparation for doing a twodimensional 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.
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 twodimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.
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.
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.
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.
Input parameters:
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.
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).
Output parameters:
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).
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 
These routines compute the threedimensional complex FFT of the
complex matrix X, and store the results in the complex matrix Y.
In FFT applications, it is customary to use zerobased subscripts; the formulas are simpler that way. So suppose the arrays are dimensioned as follows:
These routines compute the formula:
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).
The following examples are for Origin series only.
Example 1: Initialize the TABLE array in preparation for doing a threedimensional 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.
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 threedimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.
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.
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.
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.
Input parameters:
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.
Integer. Transform size in the first dimension. If n1 is not positive, the routine returns without computing a transform.
Integer. Transform size in the second dimension. If n2 is not positive, the routine returns without computing a transform.
Integer. Transform size in the third dimension. If n3 is not positive, the routine returns without computing a transform.
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.
Array of dimension (0:ldx1, 0:ldx21, 0:n31). ccfft3d
: complex array. zzfft3d
: double complex array.
Input array of values to be transformed.
Integer. The first dimension of x, as it was declared in the calling program (the leading dimension of x). ldx >= MAX(n1, 1).
Integer. The second dimension of x, as it was declared in the calling program. ldx2 >= MAX(n2, 1).
Integer. The first dimension of y, as it was declared in the calling program (the leading dimension of y). ldy >= MAX(n1, 1).
Integer. The second dimension of y, as it was declared in the calling program. ldy2 >= MAX(n2, 1).
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:
Array of dimension (0:ldy1, 0:ldy21, 0:n31). 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.
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).
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 
<tt>ccfftm/zzfftm</tt> 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:
Then column L of the output array is the FFT of column L of the input array, using the following formula for the FFT:
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 implementationdependent details.
These routines use decimationinfrequency 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.
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.
In the preceding description, it is assumed that array subscripts were zerobased, as is customary in FFT applications. Thus, the input and output arrays are declared as follows:
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:
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.
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.
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.
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.
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.
Input parameters:
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.
ccfftm
: real. zzfftm
: double precision. Each element of the output array is multiplied by scale factor after taking the Fourier transform, as defined previously. ccfftm
: real array. zzfftm
: double precision array. Input array of values to be transformed. 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 processorspecific 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), userspecified parameters are expected in isys(1) through isys(il).
Output parameters:
Array of dimension (0:ldy1, 0:lot1). 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.
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).

static 

static 

static 

static 

static 

static 

static 

static 

static 

static 
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 complextoreal transform.
It is customary in FFT applications to use zerobased subscripts; the formulas are simpler that way. For these routines, suppose that the arrays are dimensioned as follows:
Then the output array is the FFT of the input array, using the following formula for the FFT:
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.
Notice in the preceding formula that there are n real input values, and n/2 + 1 complex output values. This property is characteristic of realtocomplex 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 complextocomplex 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 realtocomplex 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 realtocomplex 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:
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 complextoreal 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.
Consider the complextoreal 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:
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 complextoreal routine, csfft
, the arrays can be dimensioned as follows:
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.
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).
In the preceding description, it is assumed that array subscripts were zerobased, as is customary in FFT applications. Thus, the input and output arrays are declared (assuming n > 0):
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:
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.
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.
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.
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.
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.
Input parameters:
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.
scfft/dzfft
returns without calculating the transform. 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. scfft
: real array of dimension (0:n1). dzfft
: double precision array of dimension (0:n1). csfft
: complex array of dimension (0:n/2). zdfft
: double complex array of dimension (0:n/2). Integer array of dimension (0:isys(0)). Use isys to specify certain processorspecific 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), userspecified parameters are expected in isys(1) through isys(il).
Output parameters:
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:n1). zdfft
: double precision array of dimension (0:n1).
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.
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.
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 

static 
<tt>scfft2d/dzfft2d</tt> computes the twodimensional 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 zerobased subscripts; the formulas are simpler that way. First the function of scfft2d
is described. Suppose the arrays are dimensioned as follows:
scfft2d
computes the formula:
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 realtocomplex transform in the first dimension, followed by the complextocomplex transform in the second dimension.
csfft2d
does the reverse. It takes the complextocomplex FFT in the second dimension, followed by the complextoreal FFT in the first dimension.
See the scfft
man page for more information about realtocomplex and complextoreal FFTs. The twodimensional analog of the conjugate formula is as follows:
Thus, you have to compute only (slightly more than) half of the output values, namely:
scfft2d
uses a routine similar to scfftm
to do a realtocomplex FFT on the columns, then uses a routine similar to ccfftm
to do a complextocomplex FFT on the rows.
csfft2d
uses a routine similar to ccfftm
to do a complextocomplex FFT on the rows, then uses a routine similar to csfftm
to do a complextoreal FFT on the columns.
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.
In the preceding description, it is assumed that array subscripts were zerobased, as is customary in FFT applications. Thus, the input and output arrays are declared:
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:
The following examples are for Origin series only.
Example 1: Initialize the TABLE array in preparation for doing a twodimensional 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.
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 twodimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.
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.
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.
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.
Input parameters:
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.
scfft2d
returns without calculating a transform. scfft2d
returns without calculating a transform. 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. Array of dimension (0:ldx1, 0:n21). scfft2d
: real array. dzfft2d
: double precision array. csfft2d
: complex array. zdfft2d
: double complex array.
Array of values to be transformed.
scfft2d, dzfft2d
: ldx >= MAX(n1, 1). csfft2d, zdfft2d
: ldx >= MAX(n1/2 + 1, 1). 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 complextoreal 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.
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 processorspecific 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), userspecified parameters are expected in isys(1) through isys(il).
Output parameters:
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.
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).
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 
These are C++ wrapper functions for the 3D realtocomplex and
complextoreal 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 threedimensional 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 zerobased subscripts; the formulas are simpler that way. First, the function of SCFFT3D
is described. Suppose the arrays are dimensioned as follows:
scfft3d
computes the formula:
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
scfft3d
is very similar in function to ccfft3d
, but it takes the realtocomplex transform in the first dimension, followed by the complextocomplex transform in the second and third dimensions.
csfft3d
does the reverse. It takes the complextocomplex FFT in the third and second dimensions, followed by the complextoreal FFT in the first dimension.
See the scfftm
man page for more information about realtocomplex and complextoreal FFTs. The three dimensional analog of the conjugate formula is as follows:
Thus, you have to compute only (slightly more than) half out the output values, namely:
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 threedimensional 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.
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 threedimensional FFT of X and store it in Y. Initialize the TABLE array, as in example 1.
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.
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.
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.
Input parameters:
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.
Integer. Transform size in the first dimension. If n1 is not positive, scfft3d
returns without computing a transform.
Integer. Transform size in the second dimension. If n2 is not positive, scfft3d
returns without computing a transform.
Integer. Transform size in the third dimension. If n3 is not positive, scfft3d
returns without computing a transform.
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.
Array of dimension (0:ldx1, 0:ldx21, 0:n31). scfft3d
: real array. dzfft3d
: double precision array. csfft3d
: complex array. zdfft3d
: double complex array.
Array of values to be transformed.
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).
Integer. The second dimension of x, as it was declared in the calling program. ldx2 >= MAX(n2, 1).
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 complextoreal 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.
Integer. The second dimension of y, as it was declared in the calling program. ldy2 >= MAX(n2, 1).
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:
Array of dimension (0:ldy1, 0:ldy21, 0:n31). 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.
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 threedimensional 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.
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 
<tt>scfftm/dzfftm</tt> 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 zerobased subscripts; the formulas are simpler that way. First, the function of scfftm
is described. Suppose that the arrays are dimensioned as follows:
Then column L of the output array is the FFT of column L of the input array, using the following formula for the FFT:
Different authors use different conventions for which transform (isign = +1 or isign = 1) is used in the realtocomplex case, and what the scale factor should be. Some adopt the convention that isign = 1 for the realtocomplex transform, and isign = 1 for the complextoreal 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 realtocomplex 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 complextoreal FFT by using isign = 1 and scale = 1.0/n.
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 realtocomplex 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 complextocomplex 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 realtocomplex 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 realtocomplex FFT routine is more efficient than the equivalent complextocomplex 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:
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 complextoreal 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.
Consider the complextoreal 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:
Thus, in the complextoreal routine, CSFFTM, the arrays can be dimensioned as follows:
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.
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.
In the preceding description, it is assumed that array subscripts were zerobased, as is customary in FFT applications. Thus, the input and output arrays are declared (for SCFFTM):
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:
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.
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.
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.
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.
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.
Input parameters:
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.
scfftm
or csfftm
returns without computing a transforms. csfftm
or scfftm
returns without computing a transforms. 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. scfftm
: real array. dzfftm
: double precision array. csfftm
: complex array. zdfftm
: double complex array. scfftm, dzfftm
: ldx >= MAX(n, 1). csfftm, zdfftm
: ldx >= MAX(n/2 + 1, 1). scfftm, dzfftm
: ldy >= MAX(n/2 + 1, 1). csfftm, zdfftm
: ldy >= MAX(n, 1). 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 processorspecific 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), userspecified parameters are expected in isys(1) through isys(il).
Output parameters:
Output array of transformed values. Dimension (0:ldy1, 0:lot1). 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.
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.

static 

static 

static 

static 

static 

static 

static 

static 