external/fftpack/fftpack5/cfft2b.F

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