Call 303-497-2400
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE CFFT1I (N, WSAVE, LENSAV, IER) INTEGER N, LENSAV, IER REAL WSAVE(LENSAV)
Send an email
Submit a help ticket
Walk-in: ML 1B Suite 55 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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE CFFT1B (N, INC, C, LENC, WSAVE, LENSAV, 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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE CFFT1F (N, INC, C, LENC, WSAVE, LENSAV, 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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE CFFT2B (LDIM, L, M, C, WSAVE, LENSAV, 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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE CFFT2F (LDIM, L, M, C, WSAVE, LENSAV, 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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE CFFTMB (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE CFFTMF (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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 (L+1)*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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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 (L+1)*M.
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE RFFTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE RFFTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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 routine
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE COSTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE COSTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE SINTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE SINTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE COSQMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE COSQMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE SINQMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * copyright (c) 2011 by UCAR * C * * C * University Corporation for Atmospheric Research * C * * C * all rights reserved * C * * C * FFTPACK version 5.1 * C * * C * A Fortran Package of Fast Fourier * C * * C * Subroutines and Example Programs * C * * C * by * C * * C * Paul Swarztrauber and Dick Valent * C * * C * of * C * * C * the National Center for Atmospheric Research * C * * C * Boulder, Colorado (80307) U.S.A. * C * * C * which is sponsored by * C * * C * the National Science Foundation * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
SUBROUTINE SINQMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 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.
FFTPACK5, Copyright (C) 2004-2011, Computational Information Systems Laboratory, University Corporation for Atmospheric Research