external/fftpack/fftpack5/zfft2f.F

   1 subroutine zfft2f ( ldim, l, m, c, wsave, lensav, work, lenwrk, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 !! ZFFT2F: complex double precision forward fast Fourier transform, 2D.
   6 !
   7 !  Discussion:
   8 !
   9 !    ZFFT2F computes the two-dimensional discrete Fourier transform of
  10 !    a complex periodic array. This transform is known as the forward
  11 !    transform or Fourier analysis, transforming from physical to
  12 !    spectral space. Routine ZFFT2F is normalized, in that a call to
  13 !    ZFFT2F followed by a call to ZFFT2B (or vice-versa) reproduces the
  14 !    original array within roundoff error.
  15 !
  16 !    On 10 May 2010, this code was modified by changing the value
  17 !    of an index into the WSAVE array.
  18 !
  19 !
  20 !
  21 !  Modified:
  22 !
  23 !    10 May 2010
  24 !
  25 !  Author:
  26 !
  27 !    Original complex single precision by Paul Swarztrauber, Richard Valent.
  28 !    Complex double precision version by John Burkardt.
  29 !
  30 !  Reference:
  31 !
  32 !    Paul Swarztrauber,
  33 !    Vectorizing the Fast Fourier Transforms,
  34 !    in Parallel Computations,
  35 !    edited by G. Rodrigue,
  36 !    Academic Press, 1982.
  37 !
  38 !    Paul Swarztrauber,
  39 !    Fast Fourier Transform Algorithms for Vector Computers,
  40 !    Parallel Computing, pages 45-63, 1984.
  41 !
  42 !  Parameters:
  43 !
  44 !    Input, integer ( kind = 4 ) LDIM, the first dimension of the array C.
  45 !
  46 !    Input, integer ( kind = 4 ) L, the number of elements to be transformed
  47 !    in the first dimension of the two-dimensional complex array C.  The value
  48 !    of L must be less than or equal to that of LDIM.  The transform is most
  49 !    efficient when L is a product of small primes.
  50 !
  51 !    Input, integer ( kind = 4 ) M, the number of elements to be transformed
  52 !    in the second dimension of the two-dimensional complex array C.  The
  53 !    transform is most efficient when M is a product of small primes.
  54 !
  55 !    Input/output, complex ( kind = 8 ) C(LDIM,M), on input, the array of two
  56 !    dimensions containing the (L,M) subarray to be transformed.  On output,
  57 !    the transformed data.
  58 !
  59 !    Input, real ( kind = 8 ) WSAVE(LENSAV). WSAVE's contents must be
  60 !    initialized with a call to ZFFT2I before the first call to routine ZFFT2F
  61 !    or ZFFT2B with transform lengths L and M.  WSAVE's contents may be re-used
  62 !    for subsequent calls to ZFFT2F and ZFFT2B having those same
  63 !    transform lengths.
  64 !
  65 !    Input, integer ( kind = 4 ) LENSAV, the dimension of the WSAVE array.
  66 !    LENSAV must be at least 2*(L+M) + INT(LOG(REAL(L)))
  67 !    + INT(LOG(REAL(M))) + 8.
  68 !
  69 !    Workspace, real ( kind = 8 ) WORK(LENWRK).
  70 !
  71 !    Input, integer ( kind = 4 ) LENWRK, the dimension of the WORK array.
  72 !    LENWRK must be at least 2*L*M.
  73 !
  74 !    Output, integer ( kind = 4 ) IER, error flag.
  75 !    0, successful exit;
  76 !    2, input parameter LENSAV not big enough;
  77 !    3, input parameter LENWRK not big enough;
  78 !    5, input parameter LDIM < L;
  79 !    20, input error returned by lower level routine.
  80 !
  81   implicit none
  82
  83   integer ( kind = 4 ) m
  84   integer ( kind = 4 ) ldim
  85   integer ( kind = 4 ) lensav
  86   integer ( kind = 4 ) lenwrk
  87
  88   complex ( kind = 8 ) c(ldim,m)
  89   integer ( kind = 4 ) ier
  90   integer ( kind = 4 ) ier1
  91   integer ( kind = 4 ) iw
  92   integer ( kind = 4 ) l
  93   real ( kind = 8 ) work(lenwrk)
  94   real ( kind = 8 ) wsave(lensav)
  95
  96   ier = 0
  97
  98   if ( ldim < l ) then
  99     ier = 5
 100     call xerfft ( 'ZFFT2F', -2 )
 101     return
 102   end if
 103
 104   if ( lensav < 2 * l + int ( log ( real ( l, kind = 8 ) ) ) + &
 105                 2 * m + int ( log ( real ( m, kind = 8 ) ) ) + 8 ) then
 106     ier = 2
 107     call xerfft ( 'ZFFT2F', 6 )
 108     return
 109   end if
 110
 111   if ( lenwrk < 2 * l * m ) then
 112     ier = 3
 113     call xerfft ( 'ZFFT2F', 8 )
 114     return
 115   end if
 116 !
 117 !  Transform the X lines of the C array.
 118 !
 119 !  On 10 May 2010, the value of IW was modified.
 120 !
 121   iw = 2 * l + int ( log ( real ( l, kind = 8 ) ) ) + 5
 122
 123   call zfftmf ( l, 1, m, ldim, c, (l-1) + ldim*(m-1) +1, wsave(iw), &
 124     2*m + int(log(real ( m, kind = 8 ))) + 4, work, 2*l*m, ier1 )
 125
 126   if ( ier1 /= 0 ) then
 127     ier = 20
 128     call xerfft ( 'ZFFT2F', -5 )
 129     return
 130   end if
 131 !
 132 !  Transform the Y lines of the C array.
 133 !
 134   iw = 1
 135
 136   call zfftmf ( m, ldim, l, 1, c, (m-1)*ldim + l, wsave(iw), &
 137     2*l + int(log(real ( l, kind = 8 ))) + 4, work, 2*m*l, ier1 )
 138
 139   if ( ier1 /= 0 ) then
 140     ier = 20
 141     call xerfft ( 'ZFFT2F', -5 )
 142     return
 143   end if
 144
 145   return
 146 end