Complex Transform Routines __________________________ CFFT1I 1D complex initialization CFFT1B 1D complex backward CFFT1F 1D complex forward CFFT2I 2D complex initialization CFFT2B 2D complex backward CFFT2F 2D complex forward CFFTMI multiple complex initialization CFFTMB multiple complex backward CFFTMF multiple complex forward Real Transform Routines _______________________ RFFT1I 1D real initialization RFFT1B 1D real backward RFFT1F 1D real forward RFFT2I 2D real initialization RFFT2B 2D real backward RFFT2F 2D real forward RFFTMI multiple real initialization RFFTMB multiple real backward RFFTMF multiple real forward Real Cosine Transform Routines ______________________________ COST1I 1D real cosine initialization COST1B 1D real cosine backward COST1F 1D real cosine forward COSTMI multiple real cosine initialization COSTMB multiple real cosine backward COSTMF multiple real cosine forward Real Sine Transform Routines ____________________________ SINT1I 1D real sine initialization SINT1B 1D real sine backward SINT1F 1D real sine forward SINTMI multiple real sine initialization SINTMB multiple real sine backward SINTMF multiple real sine forward Real Quarter-Cosine Transform Routines ______________________________________ COSQ1I 1D real quarter-cosine initialization COSQ1B 1D real quarter-cosine backward COSQ1F 1D real quarter-cosine forward COSQMI multiple real quarter-cosine initialization COSQMB multiple real quarter-cosine backward COSQMF multiple real quarter-cosine forward Real Quarter-Sine Transform Routines ____________________________________ SINQ1I 1D real quarter-sine initialization SINQ1B 1D real quarter-sine backward SINQ1F 1D real quarter-sine forward SINQMI multiple real quarter-sine initialization SINQMB multiple real quarter-sine backward SINQMF multiple real quarter-sine forward

Library FFTPACK 5.1 contains 1D, 2D, and multiple fast Fourier subroutines, written in Fortran 77, for transforming real and complex
data, real even and odd wave data, and real even and odd quarter-wave data. All of the FFTPACK 5.1 routines listed above are grouped
in triplets e.g. {CFFT1I, CFFT1F, CFFT1B}. The suffix *I* denotes initialize, *F* denotes forward (as in forward transform) and *B* denotes backward. In an application program, before calling *B* or *F* routines for the first time, or before calling them with a
different length, users must initialize an array by calling the *I* routine of the appropriate pair or triplet. Note that *I* routines
need not be called each time before a B or F routine is called.

All of the transform routines in FFTPACK 5.1 are normalized.

Error messages are written to unit 6 by routine XERFFT. The standard version of XERFFT issues an error message and halts execution, so that no FFTPACK 5.1 routine will return to the calling program with error return IER different than zero. Users may consider modifying the STOP statement in order to call system-specific exception-handling facilities.

FFTPACK 5.1 is written in standard Fortran 77 except for several instances where arrays of type REAL or COMPLEX are passed to a subroutine and used as a different type.

(1) Vectorizing the Fast Fourier Transforms, by Paul Swarztrauber, Parallel Computations, G. Rodrigue, ed., Academic Press, New York 1982.

CFFT1I 1D complex initialization

SUBROUTINE CFFT1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine CFFT1I initializes array WSAVE for use in its companion routines CFFT1B and CFFT1F. Routine CFFT1I must be called before the first call to CFFT1B or CFFT1F, and after whenever the value of integer N changes. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines CFFT1B or CFFT1F. IER = 0 successful exit = 2 input parameter LENSAV not big enough

CFFT1B 1D complex backward

SUBROUTINE CFFT1B (N, INC, C, LENC, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER N, INC, LENC, LENSAV, LENWRK, IER COMPLEX C(LENC) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFT1B computes the one-dimensional Fourier transform of a single periodic sequence within a complex array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to CFFT1B followed by a call to CFFT1F (or vice-versa) reproduces the original array subject to algorithm constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array C, of two consecutive elements within the sequence to be transformed. C Complex array of length LENC containing the sequence to be transformed. LENC Integer dimension of C array. LENC must be at least INC*(N-1) + 1. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFT1I before the first call to routine CFFT1F or CFFT1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to CFFT1F and CFFT1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*N. Output Arguments C For index J*INC+1 where J=0,...,N-1, C(J*INC+1) = N-1 SUM C(K*INC+1)*EXP(I*J*K*2*PI/N) K=0 where I=SQRT(-1). At other indices, the output value of C does not differ from input. IER = 0 successful exit = 1 input parameter LENC not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

CFFT1F 1D complex forward

SUBROUTINE CFFT1F (N, INC, C, LENC, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER N, INC, LENC, LENSAV, LENWRK, IER COMPLEX C(LENC) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFT1F computes the one-dimensional Fourier transform of a single periodic sequence within a complex array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to CFFT1F followed by a call to CFFT1B (or vice-versa) reproduces the original array subject to algorithm constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array C, of two consecutive elements within the sequence to be transformed. C Complex array of length LENC containing the sequence to be transformed. LENC Integer dimension of C array. LENC must be at least INC*(N-1) + 1.

WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFT1I before the first call to routine CFFT1F or CFFT1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to CFFT1F and CFFT1B with the same N.

LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*N. Output Arguments C For index J*INC+1 where J=0,...,N-1 (that is, for the Jth element of the sequence), C(J*INC+1) = N-1 SUM C(K*INC+1)*EXP(-I*J*K*2*PI/N) K=0 where I=SQRT(-1). At other indices, the output value of C does not differ from input. IER = 0 successful exit = 1 input parameter LENC not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

CFFT2I 2D complex initialization

SUBROUTINE CFFT2I (L, M, WSAVE, LENSAV, IER) INTEGER L, M, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 routine CFFT2I initializes real array WSAVE for use in its companion routines CFFT2F and CFFT2B for computing two- dimensional fast Fourier transforms of complex data. Prime factorizations of L and M, together with tabulations of the trigonometric functions, are computed and stored in array WSAVE. Input Arguments L Integer number of elements to be transformed in the first dimension. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension. The transform is most efficient when M is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of L and M, and also containing certain trigonometric values which will be used in routines CFFT2B or CFFT2F. WSAVE Real work array with dimension LENSAV. The WSAVE array must be initialized with a call to subroutine CFFT2I before the first call to CFFT2B or CFFT2F, and thereafter whenever the values of L, M or the contents of array WSAVE change. Using different WSAVE arrays for different transform lengths or types in the same program may reduce computation costs because the array contents can be re-used. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

CFFT2B 2D complex backward

SUBROUTINE CFFT2B (LDIM, L, M, C, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER L, M, LDIM, LENSAV, LENWRK, IER COMPLEX C(LDIM,M) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFT2B computes the two-dimensional discrete Fourier transform of a complex periodic array. This transform is known as the backward transform or Fourier synthesis, transforming from spectral to physical space. Routine CFFT2B is normalized, in that a call to CFFT2B followed by a call to CFFT2F (or vice-versa) reproduces the original array subject to algorithm constraints, roundoff error, etc. Input Arguments LDIM Integer first dimension of two-dimensional complex array C.

L Integer number of elements to be transformed in the first dimension of the two-dimensional complex array C. The value of L must be less than or equal to that of LDIM. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension of the two-dimensional complex array C. The transform is most efficient when M is a product of small primes. C Complex array of two dimensions containing the (L,M) subarray to be transformed. C's first dimension is LDIM, its second dimension must be at least M. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFT2I before the first call to routine CFFT2F or CFFT2B with transform lengths L and M. WSAVE's contents may be re-used for subsequent calls to CFFT2F and CFFT2B with the same transform lengths L and M. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. WORK Real work array. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*L*M. Output Arguments C Complex output array. For purposes of exposition, assume the index ranges of array C are defined by C(0:L-1,0:M-1). For I=0,...,L-1 and J=0,...,M-1, the C(I,J)'s are given in the traditional aliased form by L-1 M-1 C(I,J) = SUM SUM C(L1,M1)* L1=0 M1=0 EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) And in unaliased form, the C(I,J)'s are given by LF MF C(I,J) = SUM SUM C(L1,M1,K1)* L1=LS M1=MS EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) where LS= -L/2 and LF=L/2-1 if L is even; LS=-(L-1)/2 and LF=(L-1)/2 if L is odd; MS= -M/2 and MF=M/2-1 if M is even; MS=-(M-1)/2 and MF=(M-1)/2 if M is odd; and C(L1,M1) = C(L1+L,M1) if L1 is zero or negative; C(L1,M1) = C(L1,M1+M) if M1 is zero or negative; The two forms give different results when used to interpolate between elements of the sequence. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 5 input parameter L > LDIM = 20 input error returned by lower level routine

CFFT2F 2D complex forward

SUBROUTINE CFFT2F (LDIM, L, M, C, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER L, M, LDIM, LENSAV, LENWRK, IER COMPLEX C(LDIM,M) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFT2F computes the two-dimensional discrete Fourier transform of a complex periodic array. This transform is known as the forward transform or Fourier analysis, transforming from physical to spectral space. Routine CFFT2F is normalized, in that a call to CFFT2F followed by a call to CFFT2B (or vice-versa) reproduces the original array subject to algorithm constraints, roundoff error, etc. Input Arguments LDIM Integer first dimension of two-dimensional complex array C. L Integer number of elements to be transformed in the first dimension of the two-dimensional complex array C. The value of L must be less than or equal to that of LDIM. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension of the two-dimensional complex array C. The transform is most efficient when M is a product of small primes. C Complex array of two dimensions containing the (L,M) subarray to be transformed. C's first dimension is LDIM, its second dimension must be at least M.

WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFT2I before the first call to routine CFFT2F or CFFT2B with transform lengths L and M. WSAVE's contents may be re-used for subsequent calls to CFFT2F and CFFT2B having those same transform lengths.

LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8. WORK Real work array. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*L*M. Output Arguments C Complex output array. For purposes of exposition, assume the index ranges of array C are defined by C(0:L-1,0:M-1). For I=0,...,L-1 and J=0,...,M-1, the C(I,J)'s are given in the traditional aliased form by L-1 M-1 C(I,J) = 1/(L*M)*SUM SUM C(L1,M1)* L1=0 M1=0 EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) And in unaliased form, the C(I,J)'s are given by LF MF C(I,J) = 1/(L*M)*SUM SUM C(L1,M1)* L1=LS M1=MS EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) where LS= -L/2 and LF=L/2-1 if L is even; LS=-(L-1)/2 and LF=(L-1)/2 if L is odd; MS= -M/2 and MF=M/2-1 if M is even; MS=-(M-1)/2 and MF=(M-1)/2 if M is odd; and C(L1,M1) = C(L1+L,M1) if L1 is zero or negative; C(L1,M1) = C(L1,M1+M) if M1 is zero or negative; The two forms give different results when used to interpolate between elements of the sequence. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 5 input parameter L > LDIM = 20 input error returned by lower level routine

CFFTMI multiple complex initialization

SUBROUTINE CFFTMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine CFFTMI initializes array WSAVE for use in its companion routines CFFTMB and CFFTMF. Routine CFFTMI must be called before the first call to CFFTMB or CFFTMF, and after whenever the value of integer N changes. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines CFFTMB or CFFTMF. IER = 0 successful exit = 2 input parameter LENSAV not big enough

CFFTMB multiple complex backward

SUBROUTINE CFFTMB (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENC, LENSAV, LENWRK, IER COMPLEX C(LENC) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFTMB computes the one-dimensional Fourier transform of multiple periodic sequences within a complex array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to CFFTMF followed by a call to CFFTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array C. JUMP Integer increment between the locations, in array C, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array C, of two consecutive elements within the same sequence to be transformed. C Complex array containing LOT sequences, each having length N, to be transformed. C can have any number of dimensions, but the total number of locations must be at least LENC. LENC Integer dimension of C array. LENC must be at least (LOT-1)*JUMP + INC*(N-1) + 1.

WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFTMI before the first call to routine CFFTMF or CFFTMB for a given transform length N.

LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*LOT*N. Output Arguments C For index L*JUMP+J*INC+1 where J=0,...,N-1 and L=0,...,LOT-1, (that is, for the Jth element of the Lth sequence), C(L*JUMP+J*INC+1) = N-1 SUM C(L*JUMP+K*INC+1)*EXP(I*J*K*2*PI/N) K=0 where I=SQRT(-1). At other indices, the output value of C does not differ from input. IER = 0 successful exit = 1 input parameter LENC not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

CFFTMF multiple complex forward

SUBROUTINE CFFTMF (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENC, LENSAV, LENWRK, IER COMPLEX C(LENC) REAL WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine CFFTMF computes the one-dimensional Fourier transform of multiple periodic sequences within a complex array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequences from physical to spectral space. This transform is normalized since a call to CFFTMF followed by a call to CFFTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array C. JUMP Integer increment between the locations, in array C, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array C, of two consecutive elements within the same sequence to be transformed. C Complex array containing LOT sequences, each having length N, to be transformed. C can have any number of dimensions, but the total number of locations must be at least LENC. LENC Integer dimension of C array. LENC must be at least (LOT-1)*JUMP + INC*(N-1) + 1.

WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine CFFTMI before the first call to routine CFFTMF or CFFTMB for a given transform length N.

LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG(REAL(N))/LOG(2.)) + 4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least 2*LOT*N. Output Arguments C For index L*JUMP + J*INC +1 where J=0,...,N-1 and L=0,...,LOT-1, (that is, for the Jth element of the Lth sequence), C(L*JUMP+J*INC+1) = N-1 SUM C(L*JUMP+K*INC+1)*EXP(-I*J*K*2*PI/N) K=0 where I=SQRT(-1). At other indices, the output value of C does not differ from input. IER = 0 successful exit = 1 input parameter LENC not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

RFFT1I 1D real initialization

SUBROUTINE RFFT1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine RFFT1I initializes array WSAVE for use in its companion routines RFFT1B and RFFT1F. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines RFFT1B or RFFT1F. IER = 0 successful exit = 2 input parameter LENSAV not big enough

RFFT1B 1D real backward

SUBROUTINE RFFT1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFT1B computes the one-dimensional Fourier transform of a periodic sequence within a real array. This is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space.

This transform is normalized since a call to RFFT1B followed by a call to RFFT1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1) + 1.

WSAVE Real work array o length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFT1I before the first call to routine RFFT1F or RFFT1B for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). The output values of R are written over the input values. If N is even, set NH=N/2-1; then for J=0,...,N-1 R(J*INC) = R(0) + [(-1)**J*R((N-1)*INC)] NH + SUM R((2*N1-1)*INC)*COS(J*N1*2*PI/N) N1=1 NH + SUM R(2*N1*INC)*SIN(J*N1*2*PI/N) N1=1 If N is odd, set NH=(N-1)/2 and define R as above, except remove the expression in square brackets []. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough

RFFT1F 1D real forward

SUBROUTINE RFFT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFT1F computes the one-dimensional Fourier transform of a periodic sequence within a real array. This is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space.

This transform is normalized since a call to RFFT1F followed by a call to RFFT1B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1) + 1.

WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFT1I before the first call to routine RFFT1F or RFFT1B for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). Then N-1 R(0) = SUM R(N1*INC)/N N1=0 If N is even, set NH=N/2-1; if N is odd set NH=(N-1)/2; then for J=1,...,NH R((2*J-1)*INC) = N-1 2.*SUM (R(N1*INC)*COS(J*N1*2*PI/N)/N N1=0 and R(2*J*INC) = N-1 2.*SUM (R(N1*INC)*SIN(J*N1*2*PI/N)/N N1=0 Also if N is even then R((N-1)*INC) = N-1 SUM (-1)**N1*R(N1*INC)/N N1=0 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough

RFFT2I 2D real initialization

SUBROUTINE RFFT2I (L, M, WSAVE, LENSAV, IER) INTEGER L, M, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 routine RFFT2I initializes real array WSAVE for use in its companion routines RFFT2F and RFFT2B for computing the two- dimensional fast Fourier transform of real data. Prime factorizations of L and M, together with tabulations of the trigonometric functions, are computed and stored in array WSAVE. RFFT2I must be called prior to the first call to RFFT2F or RFFT2B. Separate WSAVE arrays are required for different values of L or M. Input Arguments L Integer number of elements to be transformed in the first dimension. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension. The transform is most efficient when M is a product of small primes. LENSAV Integer number of elements in the WSAVE array. LENSAV must be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) + 2*INT(LOG(REAL(M))/LOG(2.)) +12. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of L and M, and also containing certain trigonometric values which will be used in routines RFFT2B or RFFT2F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

RFFT2B 2D real backward

SUBROUTINE RFFT2B (LDIM, L, M, R, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LDIM, L, M, LENSAV, LENWRK, IER REAL R(LDIM,M), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFT2B computes the two-dimensional discrete Fourier transform of the complex Fourier coefficients a real periodic array. This transform is known as the backward transform or Fourier synthesis, transforming from spectral to physical space. Routine RFFT2B is normalized: a call to RFFT2B followed by a call to RFFT2F (or vice-versa) reproduces the original array subject to algorithmic sonstraints, roundoff error, etc. Input Arguments LDIM Integer first dimension of the two-dimensional real array R, which must be at least L. L Integer number of elements to be transformed in the first dimension of the two-dimensional real array R. The value of L must be less than or equal to that of LDIM. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension of the two-dimensional real array R. The transform is most efficient and accurate when M is a product of small primes. R A real L by M array containing the spectral coefficients of a real L by M array that are stored as described in the documentation of subroutine rfft2f. THe first dimension is LDIM which must be at least L. The second dimension must be at least M. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFT2I before the first call to routine RFFT2F or RFFT2B. WSAVE's contents may be re-used for subsequent calls to RFFT2F and RFFT2B with the same L and M. LENSAV Integer number of elements in the WSAVE array. LENSAV must be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) + 2*INT(LOG(REAL(M))/LOG(2.)) +12. WORK Real array of dimension LENWRK, where LENWRK is defined below. WORK provides workspace, and its contents need not be saved between calls to routines RFFT2B and RFFT2F. LENWRK Integer number of elements in the WORK array. LENWRK must be at least LDIM*M. Output Arguments R A real L by M array. If the full transform c is reconstructed using subroutine r2c(ldim,lcdim,l,m,r,c) then for i=0,...,l-1 and j=0,...,m-1 L-1 M-1 R(I,J) = SUM SUM C(L1,M1) L1=0 M1=0 *EXP(SQRT(-1)*2*PI*(I*L1/L+J*M1/M)) or using the conjugate symmetry c(i,j)=c(l-1,m-j) this can be written in terms of c(i,j), i=0,...,(L+1)/2 and j=0,...,m as: M-1 R(I,J) = REAL[ SUM C(0,M1)*EXP(SQRT(-1)*2*PI*J*M1/M ] M1=0 (L+1)/2-1 M-1 + 2*REAL[ SUM SUM C(L1,M1)* L1=1 M1=0 *EXP(SQRT(-1)*2*PI*(I*L1/L+J*M1/M)) ] If L is even then add M-1 + REAL[ SUM (-1)**I*C(L/2,M1)*EXP(SQRT(-1)*2*PI*J*M1/M) ] M1=0 c(i,j) = a(i,j)+i*b(I,J) are contained in the real output array r(i,j) except for c(0,m-j) and c(l,m-j) which can be obtained from c(0,m-j) = c(0,j) and c(l,j) = c(l.m-j) = c(l,j) IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 6 input parameter LDIM is less than 2*INT((L+1)/2) = 20 input error returned by lower level routine

RFFT2F 2D real forward

SUBROUTINE RFFT2F (LDIM, L, M, R, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER LDIM, L, M, LENSAV, LENWRK, IER REAL R(LDIM,M), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine RFFT2F computes the two-dimensional discrete Fourier transform of a real periodic array. This transform is known as the forward transform or Fourier analysis, transforming from physical to spectral space. Routine RFFT2F is normalized: a call to RFFT2F followed by a call to RFFT2B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LDIM Integer first dimension of the two-dimensional real array R, which must be at least L. L Integer number of elements to be transformed in the first dimension of the two-dimensional real array R. The value of L must be less than or equal to that of LDIM. The transform is most efficient when L is a product of small primes. M Integer number of elements to be transformed in the second dimension of the two-dimensional real array R. The transform is most efficient when M is a product of small primes. R Real array of two dimensions containing the L-by-M subarray to be transformed. R's first dimension is LDIM and its second dimension must be at least as large as M. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFT2I before the first call to routine RFFT2F or RFFT2B. WSAVE's contents may be re-used for subsequent calls to RFFT2F and RFFT2B as long as L and M remain unchanged. LENSAV Integer number of elements in the WSAVE array. LENSAV must be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) + 2*INT(LOG(REAL(M))/LOG(2.)) +12. WORK Real array of dimension LENWRK which is defined below. WORK provides workspace, and its contents need not be saved between calls to routines RFFT2F and RFFT2B. LENWRK Integer number of elements in the WORK array. LENWRK must be at least M*(L+1). Output Arguments R Real output array of two dimensions. The full complex transform of r(i,j) is given by: L-1 M-1 C(I,J) = 1/(L*M)*SUM SUM R(L1,M1)* L1=0 M1=0 EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M)) The complex transform of a real array has conjugate symmetry. That is, c(i,j) = congugate c(l-i,m-j) so only half the transform is computed and packed back into the original array R. Examples: Let a(i,j) = re[c(i,j)] and b(i,j) = Im[c(i,j)] then following the forward transform For l=m=4 a(0,0) a(0,1) b(0,1) a(0,2) r(i,j) = a(1,0) a(1,1) a(1,2) a(1,3) b(1,0) b(1,1) b(1,2) b(1,3) a(2,0) a(2,1) b(2,1) a(2,2) For l=m=5 a(0,0) a(0,1) b(0,1) a(0,2),b(0,2) a(1,0) a(1,1) a(1,2) a(1,3),a(1,4) r(i,j) = b(1,0) b(1,1) b(1,2) b(1,3),b(1,4) a(2,0) a(2,1) a(2,2) a(2,3),a(2,4) b(2,0) b(2,1) b(2,2) b(2,3),b(2,4) The remaining c(i,j) for i=int((L+1)/2)+1,..,L and m=0,...m-1 can be obtained via the conjugate symmetry, which also implies that c(0,j) = conjugate c(0,m-j) and for even l, c(l/2,0) = conjugate c(l/2,m-j). The full complex transform c(i,j), i=1,...,L and j=1,...,M can also be extracted using subroutine r2c(ldim,lcdim,l,m,r,c) where lcdim is the first dimension of the complex array c, which must be greater than or equal to l. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 6 input parameter LDIM is less than 2*INT((L+1)/2) = 20 input error returned by lower level routine

RFFTMI multiple real initialization

SUBROUTINE RFFTMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine RFFTMI initializes array WSAVE for use in its companion routines RFFTMB and RFFTMF. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines RFFTMB or RFFTMF. IER = 0 successful exit = 2 input parameter LENSAV not big enough

RFFTMB multiple real backward

SUBROUTINE RFFTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV) ,WORK(LENWRK)

FFTPACK 5.1 routine RFFTMB computes the one-dimensional Fourier transform of multiple periodic sequences within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to RFFTMB followed by a call to RFFTMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1) + 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFTMI before the first call to routine RFFTMF or RFFTMB for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. If N is even, set NH=N/2-1; then for I=0,...,LOT-1 and J=0,...,N-1 R(I*JUMP+J*INC) = R(I*JUMP) + [(-1)**J*R(I*JUMP+(N-1)*INC)] NH + SUM R(I*JUMP+(2*N1-1)*INC)*COS(J*N1*2*PI/N) N1=1 NH + SUM R(I*JUMP+2*N1*INC)*SIN(J*N1*2*PI/N) N1=1 If N is odd, set NH=(N-1)/2 and define R as above, except remove the expression in square brackets []. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

RFFTMF multiple real forward

SUBROUTINE RFFTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV) ,WORK(LENWRK)

FFTPACK 5.1 routine RFFTMF computes the one-dimensional Fourier transform of multiple periodic sequences within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequences from physical to spectral space. This transform is normalized since a call to RFFTMF followed by a call to RFFTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1) + 1. WSAVE Real work array o length LENSAV. WSAVE's contents must be initialized with a call to subroutine RFFTMI before the first call to routine RFFTMF or RFFTMB for a given transform length N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). Then for I=0,...,LOT-1 N-1 R(I*JUMP) = SUM R(I*JUMP+N1*INC)/N N1=0 If N is even, set NH=N/2-1; if N is odd set NH=(N-1)/2; then for J=1,...,NH R(I*JUMP+(2*J-1)*INC) = N-1 2.*SUM (R(I*JUMP+N1*INC)*COS(J*N1*2*PI/N)/N N1=0 and R(I*JUMP+2*J*INC) = N-1 2.*SUM (R(I*JUMP+N1*INC)*SIN(J*N1*2*PI/N)/N N1=0 Also if N is even then R(I*JUMP+(N-1)*INC) = N-1 SUM (-1)**N1*R(I*JUMP+N1*INC)/N N1=0 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

COST1I 1D real cosine initialization

SUBROUTINE COST1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine COST1I initializes array WSAVE for use in its companion routines COST1F and COST1B. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines COST1B or COST1F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

COST1B 1D real cosine backward

SUBROUTINE COST1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COST1B computes the one-dimensional Fourier transform of an even sequence within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to COST1B followed by a call to COST1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COST1I before the first call to routine COST1F or COST1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COST1F and COST1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N-1. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). The output values of R are written over the input values. For J=0,...,N-1 R(J*INC) = N-1 SUM R(N1*INC)*COS(J*N1*PI/(N-1)) N1=0 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

COST1F 1D real cosine forward

SUBROUTINE COST1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COST1F computes the one-dimensional Fourier transform of an even sequence within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to COST1F followed by a call to COST1B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. INC Integer increment between the locations, in array R of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COST1I before the first call to routine COST1F or COST1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COST1F and COST1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N-1. Output Arguments R Real output array R. For purposes of this exposition, assume R's range of indices is given by R(0:(N-1)*INC) and that the input values of R are denoted by X. The output values of R are written over the input values X as follows: CASE N=1 -------- R(0) = X(0) CASE N>1 -------- For J=0 R(0) = 0.5*X(0)/(N-1) N-2 + SUM R(N1*INC)/(N-1) N1=1 + 0.5*X((N-1)*INC)/(N-1) For J=1,...,N-2 R(J*INC) = R(0)/(N-1) N-2 + SUM 2.0*(X(N1*INC)*COS(J*N1*PI/(N-1)))/(N-1) N1=1 + ((-1)**J)*X((N-1)*INC)/(N-1) R((N-1)*INC) = 0.5*X(0)/(N-1) N-2 + SUM R(N1*INC)*((-1)**N1)/(N-1) N1=1 + 0.5*((-1)**(N-1))*X((N-1)*INC)/(N-1) IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routineCOSTMI multiple real cosine initialization

SUBROUTINE COSTMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine COSTMI initializes array WSAVE for use in its companion routines COSTMF and COSTMB. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4 Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines COSTMB or COSTMF. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

COSTMB multiple real cosine backward

SUBROUTINE COSTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSTMB computes the one-dimensional Fourier transform of multiple even sequences within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to COSTMB followed by a call to COSTMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COSTMI before the first call to routine COSTMF or COSTMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSTMF and COSTMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*(N+1). Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=0,...,N-1 R(I*JUMP+J*INC) = N-1 SUM R(I*JUMP+N1*INC)*COS(J*N1*PI/(N-1)) N1=0 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent, otherwise at least one array element mistakenly is transformed more than once.

COSTMF multiple real cosine forward

SUBROUTINE COSTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSTMF computes the one-dimensional Fourier transform of multiple even sequences within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequences from physical to spectral space. This transform is normalized since a call to COSTMF followed by a call to COSTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N-1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COSTMI before the first call to routine COSTMF or COSTMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSTMF and COSTMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*(N+1). Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 R(I*JUMP) = 0.5*X(I*JUMP)/(N-1) N-2 + SUM R(I*JUMP+*N1*INC)/(N-1) N1=1 + 0.5*X(I*JUMP+(N-1)*INC)/(N-1) For I=0,...,LOT-1 and J=1,...,N-2 R(I*JUMP+J*INC) = R(I*JUMP)/(N-1) N-2 + SUM 2.0*(X(I*JUMP+*N1*INC)*COS(J*N1*PI/(N-1)))/(N-1) N1=1 + ((-1)**J)*X(I*JUMP+(N-1)*INC)/(N-1) For I=0,...,LOT-1 R(I*JUMP+(N-1)*INC) = 0.5*X(I*JUMP)/(N-1) N-2 + SUM R(I*JUMP+*N1*INC)*((-1)**N1)/(N-1) N1=1 + 0.5*((-1)**(N-1))*X(I*JUMP+(N-1)*INC)/(N-1) IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent, otherwise at least one array element mistakenly is transformed more than once.

SINT1I 1D real sine initialization

SUBROUTINE SINT1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine SINT1I initializes array WSAVE for use in its companion routines SINT1F and SINT1B. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines SINT1B or SINT1F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SINT1B 1D real sine backward

SUBROUTINE SINT1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINT1B computes the one-dimensional Fourier transform of an odd sequence within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to SINT1B followed by a call to SINT1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINT1I before the first call to routine SINT1F or SINT1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINT1F and SINT1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. Must be at least 2*N+2. Output Arguments R Real output array. For purposes of exposition, assume R's range of indices is given by R(INC:N*INC). The output values of R are written over the input values. For J=1,...,N R(J*INC) = N SUM R(N1*INC)*SIN(J*N1*PI/(N+1)) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SINT1F 1D real sine forward

SUBROUTINE SINT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINT1F computes the one-dimensional Fourier transform of an odd sequence within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to SINT1F followed by a call to SINT1B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINT1I before the first call to routine SINT1F or SINT1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINT1F and SINT1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. Must be at least 2*N+2. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:(N-1)*INC). The output values of R are written over the input values. For J=1,...,N R(J*INC) = N SUM 2.*R(N1*INC)*SIN(J*N1*PI/(N+1))/(N+1) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SINTMI multiple real sine initialization

SUBROUTINE SINTMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine SINTMI initializes array WSAVE for use in its companion routines SINTMF and SINTMB. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines SINTMB or SINTMF. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SINTMB multiple real sine backward

SUBROUTINE SINTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINTMB computes the one-dimensional Fourier transform of multiple odd sequences within a real array. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to SINTMB followed by a call to SINTMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences. N Integer length of each sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINTMI before the first call to routine SINTMF or SINTMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINTMF and SINTMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*(2*N+4). Output Arguments R Real output array. For purposes of exposition, assume R's range of indices is given by R(INC:(LOT-1)*JUMP+N*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=1,...,N R(I*JUMP+J*INC) = N SUM R(I*JUMP+*N1*INC)*SIN(J*N1*PI/(N+1)) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SINTMF multiple real sine forward

SUBROUTINE SINTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINTMF computes the one-dimensional Fourier transform of multiple odd sequences within a real array. This transform is referred to as the forward transform or Fourier analysis, transforming the sequences from physical to spectral space. This transform is normalized since a call to SINTMF followed by a call to SINTMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N+1 is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINTMI before the first call to routine SINTMF or SINTMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINTMF and SINTMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*(2*N+4). Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=1,...,N R(I*JUMP+J*INC) = N SUM 2.*R(I*JUMP+*N1*INC)*SIN(J*N1*PI/(N+1))/(N+1) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

COSQ1I 1D real quarter-cosine initialization

SUBROUTINE COSQ1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine COSQ1I initializes array WSAVE for use in its companion routines COSQ1F and COSQ1B. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines COSQ1B or COSQ1F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

COSQ1B 1D real quarter-cosine backward

SUBROUTINE COSQ1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSQ1B computes the one-dimensional Fourier transform of a sequence which is a cosine series with odd wave numbers. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to COSQ1B followed by a call to COSQ1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer number of elements to be transformed in the sequence. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine COSQ1I before the first call to routine COSQ1F or COSQ1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSQ1F and COSQ1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N. Output Arguments R Real output array. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). The output values of R are written over the input values. For J=0,...,N-1 R(J*INC) = N-1 SUM R(N1*INC)*COS(J*(2*N1+1)*PI/(2*N)) N1=0 WSAVE Contains values initialized by subroutine COSQ1I that must not be destroyed between calls to routine COSQ1F or COSQ1B. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

COSQ1F 1D real quarter-cosine forward

SUBROUTINE COSQ1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSQ1F computes the one-dimensional Fourier transform of a sequence which is a cosine series with odd wave numbers. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to COSQ1F followed by a call to COSQ1B (or vice-versa) reproduces the original array subject to algorithmic constrain, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine COSQ1I before the first call to routine COSQ1F or COSQ1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSQ1F and COSQ1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(N-1)*INC). The output values of R are written over the input values. For J=0,...,N-1 R(J*INC) = R(0)/N N-1 + SUM 2.*R(N1*INC)*COS((2*J+1)*N1*PI/(2*N))/N N1=1 WSAVE Contains values initialized by subroutine COSQ1I that must not be destroyed between calls to routine COSQ1F or COSQ1B. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

COSQMI multiple real quarter-cosine initialization

SUBROUTINE COSQMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine COSQMI initializes array WSAVE for use in its companion routines COSQMF and COSQMB. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines COSQMB or COSQMF. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

COSQMB multiple real quarter-cosine backward

SUBROUTINE COSQMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSQMB computes the one-dimensional Fourier transform of multiple sequences, each of which is a cosine series with odd wave numbers. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to COSQMB followed by a call to COSQMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array with dimension LENSAV. WSAVE's contents must be initialized with a call to subroutine COSQMI before the first call to routine COSQMF or COSQMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSQMF and COSQMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=0,...,N-1 R(I*JUMP+J*INC) = N-1 SUM R(I*JUMP+N1*INC)*COS(J*(2*N1+1)*PI/(2*N)) N1=0 WSAVE Contains values initialized by subroutine COSQMI that must not be destroyed between calls to routine COSQMF or COSQMB. IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent, otherwise at least one array element mistakenly is transformed more than once.

COSQMF multiple real quarter-cosine forward

SUBROUTINE COSQMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine COSQMF computes the one-dimensional Fourier transform of multiple sequences within a real array, where each of the sequences is a cosine series with odd wave numbers. This transform is referred to as the forward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to COSQMF followed by a call to COSQMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array o length LENSAV. WSAVE's contents must be initialized with a call to subroutine COSQMI before the first call to routine COSQMF or COSQMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to COSQMF and COSQMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real array of dimension LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(0:(LOT-1)*JUMP+(N-1)*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=0,...,N-1 R(I*JUMP+J*INC) = R(I*JUMP)/N N-1 + SUM 2.*R(I*JUMP+*N1*INC)*COS((2*J+1)*N1*PI/(2*N))/N N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent, otherwise at least one array element mistakenly is transformed more than once.

SINQ1I 1D real quarter-sine initialization

SUBROUTINE SINQ1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine SINQ1I initializes array WSAVE for use in its companion routines SINQ1F and SINQ1B. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines SINQ1B or SINQ1F. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SINQ1B 1D real quarter-sine backward

SUBROUTINE SINQ1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINQ1B computes the one-dimensional Fourier transform of a sequence which is a sine series with odd wave numbers. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequence from spectral to physical space. This transform is normalized since a call to SINQ1B followed by a call to SINQ1F (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINQ1I before the first call to routine SINQ1F or SINQ1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINQ1F and SINQ1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:N*INC). The output values of R are written over the input values. For J=1,...,N R(J*INC) = N SUM R(N1*INC)*SIN(J*(2*N1-1)*PI/(2*N)) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SINQ1F 1D real quarter-sine forward

SUBROUTINE SINQ1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER) INTEGER N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINQ1F computes the one-dimensional Fourier transform of a sequence which is a sine series of odd wave numbers. This transform is referred to as the forward transform or Fourier analysis, transforming the sequence from physical to spectral space. This transform is normalized since a call to SINQ1F followed by a call to SINQ1B (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments N Integer length of the sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the sequence. R Real array of length LENR containing the sequence to be transformed. LENR Integer dimension of R array. LENR must be at least INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINQ1I before the first call to routine SINQ1F or SINQ1B for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINQ1F and SINQ1B with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:N*INC). The output values of R are written over the input values. For J=1,...,N R(J*INC) = N-1 + SUM (2.*R(N1*INC)*SIN(((2*J-1)*N1*PI/(2*N)))/N N1=1 + ((-1)**(J+1))*R(N*INC)/N IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 20 input error returned by lower level routine

SINQMI multiple real quarter-sine initialization

SUBROUTINE SINQMI (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)

FFTPACK 5.1 subroutine SINQMI initializes array WSAVE for use in its companion routines SINQMF and SINQMB. The prime factor- ization of N together with a tabulation of the trigonometric functions are computed and stored in array WSAVE. Separate WSAVE arrays are required for different values of N. Input Arguments N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. Output Arguments WSAVE Real work array with dimension LENSAV, containing the prime factors of N and also containing certain trigonometric values which will be used in routines SINQMB or SINQMF. IER Integer error return = 0 successful exit = 2 input parameter LENSAV not big enough = 20 input error returned by lower level routine

SINQMB multiple real quarter-sine backward

SUBROUTINE SINQMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINQMB computes the one-dimensional Fourier transform of multiple sequences within a real array, where each of the sequences is a sine series with odd wave numbers. This transform is referred to as the backward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to SINQMB followed by a call to SINQMF (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINQMI before the first call to routine SINQMF or SINQMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINQMF and SINQMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:(LOT-1)*JUMP+N*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=1,...,N R(I*JUMP+J*INC) = N SUM R(I*JUMP+N1*INC)*SIN(J*(2*N1-1)*PI/(2*N)) N1=1 IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.

SINQMF multiple real quarter-sine forward

SUBROUTINE SINQMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 1 WORK, LENWRK, IER) INTEGER LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER REAL R(LENR), WSAVE(LENSAV), WORK(LENWRK)

FFTPACK 5.1 routine SINQMF computes the one-dimensional Fourier transform of multiple sequences within a real array, where each sequence is a sine series with odd wave numbers. This transform is referred to as the forward transform or Fourier synthesis, transforming the sequences from spectral to physical space. This transform is normalized since a call to SINQMF followed by a call to SINQMB (or vice-versa) reproduces the original array subject to algorithmic constraints, roundoff error, etc. Input Arguments LOT Integer number of sequences to be transformed within array R. JUMP Integer increment between the locations, in array R, of the first elements of two consecutive sequences to be transformed. N Integer length of each sequence to be transformed. The transform is most efficient when N is a product of small primes. INC Integer increment between the locations, in array R, of two consecutive elements within the same sequence. R Real array containing LOT sequences, each having length N. R can have any number of dimensions, but the total number of locations must be at least LENR. LENR Integer dimension of R array. LENR must be at least (LOT-1)*JUMP + INC*(N-1)+ 1. WSAVE Real work array of length LENSAV. WSAVE's contents must be initialized with a call to subroutine SINQMI before the first call to routine SINQMF or SINQMB for a given transform length N. WSAVE's contents may be re-used for subsequent calls to SINQMF and SINQMB with the same N. LENSAV Integer dimension of WSAVE array. LENSAV must be at least 2*N + INT(LOG (REAL(N))/LOG(2.)) +4. WORK Real work array of dimension at least LENWRK. LENWRK Integer dimension of WORK array. LENWRK must be at least LOT*N. Output Arguments R Real output array R. For purposes of exposition, assume R's range of indices is given by R(INC:(LOT-1)*JUMP+N*INC). The output values of R are written over the input values. For I=0,...,LOT-1 and J=1,...,N R(I*JUMP+J*INC) = N-1 + SUM (2.*R(I*JUMP+*N1*INC)*SIN(((2*J-1)*N1*PI/(2*N)))/N N1=1 + ((-1)**(J+1))*R(I*JUMP+N*INC)/N IER Integer error return = 0 successful exit = 1 input parameter LENR not big enough = 2 input parameter LENSAV not big enough = 3 input parameter LENWRK not big enough = 4 input parameters INC,JUMP,N,LOT are not consistent. = 20 input error returned by lower level routine The parameters integers INC, JUMP, N and LOT are consistent if equality I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N and J1,J2 < LOT implies I1=I2 and J1=J2. For multiple FFTs to execute correctly, input variables INC, JUMP, N and LOT must be consistent ... otherwise at least one array element mistakenly is transformed more than once.