external/fftpack/fftpack5/rfft2f.F

   1 subroutine rfft2f ( ldim, l, m, r, wsave, lensav, work, lenwrk, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 ! RFFT2F: real single precision forward fast Fourier transform, 2D.
   6 !
   7 !  Discussion:
   8 !
   9 !    RFFT2F computes the two-dimensional discrete Fourier transform of a
  10 !    real periodic array.  This transform is known as the forward transform
  11 !    or Fourier analysis, transforming from physical to spectral space.
  12 !    Routine RFFT2F is normalized: a call to RFFT2F followed by a call to
  13 !    RFFT2B (or vice-versa) reproduces the original array within roundoff
  14 !    error.
  15 !
  16 !
  17 !    Copyright (C) 1995-2004, Scientific Computing Division,
  18 !    University Corporation for Atmospheric Research
  19 !
  20 !  Modified:
  21 !
  22 !    26 March 2005
  23 !
  24 !  Author:
  25 !
  26 !    Paul Swarztrauber
  27 !    Richard Valent
  28 !
  29 !  Reference:
  30 !
  31 !    Paul Swarztrauber,
  32 !    Vectorizing the Fast Fourier Transforms,
  33 !    in Parallel Computations,
  34 !    edited by G. Rodrigue,
  35 !    Academic Press, 1982.
  36 !
  37 !    Paul Swarztrauber,
  38 !    Fast Fourier Transform Algorithms for Vector Computers,
  39 !    Parallel Computing, pages 45-63, 1984.
  40 !
  41 !  Parameters:
  42 !
  43 !    Input, integer ( kind = 4 ) LDIM, the first dimension of the 2D real
  44 !    array R, which must be at least 2*(L/2+1).
  45 !
  46 !    Input, integer ( kind = 4 ) L, the number of elements to be transformed
  47 !    in the first dimension of the two-dimensional real array R.  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 real array R.  The
  53 !    transform is most efficient when M is a product of small primes.
  54 !
  55 !    Input/output, real ( kind = 4 ) R(LDIM,M), the real array of two
  56 !    dimensions.  On input, containing the L-by-M physical data to be
  57 !    transformed.  On output, the spectral coefficients.
  58 !
  59 !    Input, real ( kind = 4 ) WSAVE(LENSAV).  WSAVE's contents must be
  60 !    initialized with a call to RFFT2I before the first call to routine RFFT2F
  61 !    or RFFT2B with lengths L and M.  WSAVE's contents may be re-used for
  62 !    subsequent calls to RFFT2F and RFFT2B with the same transform lengths.
  63 !
  64 !    Input, integer ( kind = 4 ) LENSAV, the number of elements in the WSAVE
  65 !    array.  LENSAV must be at least L + M + INT(LOG(REAL(L)))
  66 !    + INT(LOG(REAL(M))) + 8.
  67 !
  68 !    Workspace, real ( kind = 4 ) WORK(LENWRK), provides workspace, and its
  69 !    contents need not be saved between calls to routines RFFT2F and RFFT2B.
  70 !
  71 !    Input, integer ( kind = 4 ) LENWRK, the number of elements in the WORK
  72 !    array.  LENWRK must be at least LDIM*M.
  73 !
  74 !    Output, integer ( kind = 4 ) IER, the error flag.
  75 !    0, successful exit;
  76 !    2, input parameter LENSAV not big enough;
  77 !    3, input parameter LENWRK not big enough;
  78 !    6, input parameter LDIM < 2*(L+1);
  79 !    20, input error returned by lower level routine.
  80 !
  81   implicit none
  82
  83   integer ( kind = 4 ) ldim
  84   integer ( kind = 4 ) lensav
  85   integer ( kind = 4 ) lenwrk
  86   integer ( kind = 4 ) m
  87
  88   integer ( kind = 4 ) i
  89   integer ( kind = 4 ) ier
  90   integer ( kind = 4 ) ier1
  91   integer ( kind = 4 ) j
  92   integer ( kind = 4 ) l
  93   integer ( kind = 4 ) lwsav
  94   integer ( kind = 4 ) mwsav
  95   real ( kind = 4 ) work(lenwrk)
  96   real ( kind = 4 ) wsave(lensav)
  97   real ( kind = 4 ) r(ldim,m)
  98
  99   ier = 0
 100 !
 101 !  Verify LENSAV.
 102 !
 103   lwsav = l + int ( log ( real ( l, kind = 4 ) ) ) + 4
 104   mwsav = 2 * m + int ( log ( real ( m, kind = 4 ) ) ) + 4
 105
 106   if ( lensav < lwsav + mwsav ) then
 107     ier = 2
 108     call xerfft ( 'rfft2f', 6 )
 109     return
 110   end if
 111 !
 112 !  Verify LENWRK.
 113 !
 114   if ( lenwrk < 2 * ( l / 2 + 1 ) * m ) then
 115     ier = 3
 116     call xerfft ( 'rfft2f', 8 )
 117     return
 118   end if
 119 !
 120 !  Verify LDIM is as big as L.
 121 !
 122   if ( ldim < 2 * ( l / 2 + 1 ) ) then
 123     ier = 5
 124     call xerfft ( 'rfft2f', -6 )
 125     return
 126   end if
 127 !
 128 !  Transform first dimension of array.
 129 !
 130   call rfftmf ( m, ldim, l, 1, r, m*ldim, wsave(1), &
 131     l+int(log( real ( l, kind = 4 )))+4, work,2*(l/2+1)*m, ier1 )
 132
 133   if ( ier1 /= 0 ) then
 134      ier = 20
 135      call xerfft ( 'rfft2f', -5 )
 136      return
 137   end if
 138 !
 139 !  Reshuffle to add in Nyquist imaginary components.
 140 !
 141   do j = 1, m
 142     if ( mod ( l, 2 ) == 0 ) then
 143       r(l+2,j) = 0.0E+00
 144     end if
 145     do i = l, 2, -1
 146       r(i+1,j) = r(i,j)
 147     end do
 148     r(2,j) = 0.0E+00
 149   end do
 150 !
 151 !  Transform second dimension of array.
 152 !
 153   call cfftmf ( l/2+1, 1, m, ldim/2, r, m*ldim/2, &
 154     wsave(l+int(log( real ( l, kind = 4 )))+5), &
 155     2*m+int(log( real ( m, kind = 4 )))+4, work, 2*(l/2+1)*m, ier1 )
 156
 157   if ( ier1 /= 0 ) then
 158     ier = 20
 159     call xerfft ( 'rfft2f', -5 )
 160     return
 161   end if
 162
 163   return
 164 end