external/fftpack/fftpack5/rfft2b.F

   1 subroutine rfft2b ( ldim, l, m, r, wsave, lensav, work, lenwrk, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 !! RFFT2B: real single precision backward fast Fourier transform, 2D.
   6 !
   7 !  Discussion:
   8 !
   9 !    RFFT2B computes the two-dimensional discrete Fourier transform of the
  10 !    complex Fourier coefficients a real periodic array.  This transform is
  11 !    known as the backward transform or Fourier synthesis, transforming from
  12 !    spectral to physical space.  Routine RFFT2B is normalized: a call to
  13 !    RFFT2B followed by a call to RFFT2F (or vice-versa) reproduces the
  14 !    original array within roundoff 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 of
  48 !    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 transform
  53 !    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, R contains the L/2+1-by-M complex subarray of
  57 !    spectral coefficients, on output, the physical 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 !    L and M.
  64 !
  65 !    Input, integer ( kind = 4 ) LENSAV, the number of elements in the WSAVE
  66 !    array.  LENSAV must be at least L + M + INT(LOG(REAL(L)))
  67 !    + INT(LOG(REAL(M))) + 8.
  68 !
  69 !    Workspace, real ( kind = 4 ) WORK(LENWRK).  WORK provides workspace, and
  70 !    its contents need not be saved between calls to routines RFFT2B and RFFT2F.
  71 !
  72 !    Input, integer ( kind = 4 )  LENWRK, the number of elements in the WORK
  73 !    array.  LENWRK must be at least LDIM*M.
  74 !
  75 !    Output, integer ( kind = 4 ) IER, the error flag.
  76 !    0, successful exit;
  77 !    2, input parameter LENSAV not big enough;
  78 !    3, input parameter LENWRK not big enough;
  79 !    6, input parameter LDIM < 2*(L/2+1);
  80 !    20, input error returned by lower level routine.
  81 !
  82   implicit none
  83
  84   integer ( kind = 4 ) ldim
  85   integer ( kind = 4 ) lensav
  86   integer ( kind = 4 ) lenwrk
  87   integer ( kind = 4 ) m
  88
  89   integer ( kind = 4 ) i
  90   integer ( kind = 4 ) ier
  91   integer ( kind = 4 ) ier1
  92   integer ( kind = 4 ) j
  93   integer ( kind = 4 ) l
  94   integer ( kind = 4 ) lwsav
  95   integer ( kind = 4 ) mwsav
  96   real ( kind = 4 ) work(lenwrk)
  97   real ( kind = 4 ) wsave(lensav)
  98   real ( kind = 4 ) r(ldim,m)
  99
 100   ier = 0
 101 !
 102 !  Verify LENSAV.
 103 !
 104   lwsav = l + int ( log ( real ( l, kind = 4 ) ) ) + 4
 105   mwsav = 2 * m + int ( log ( real ( m, kind = 4 ) ) ) + 4
 106
 107   if ( lensav < lwsav + mwsav ) then
 108     ier = 2
 109     call xerfft ( 'rfft2b', 6 )
 110     return
 111   end if
 112 !
 113 !  Verify LENWRK.
 114 !
 115   if ( lenwrk < 2 * ( l / 2 + 1 ) * m ) then
 116     ier = 3
 117     call xerfft ( 'rfft2b', 8 )
 118     return
 119   end if
 120 !
 121 !  Verify LDIM is as big as L.
 122 !
 123   if ( ldim < 2 * ( l / 2 + 1 ) ) then
 124     ier = 5
 125     call xerfft ( 'rfft2b', -6 )
 126     return
 127   end if
 128 !
 129 !  Transform second dimension of array.
 130 !
 131   call cfftmb ( l/2+1, 1, m, ldim/2, r, m*ldim/2, &
 132     wsave(l+int(log( real ( l, kind = 4 )))+5), &
 133     2*m+int(log( real ( m, kind = 4 )))+4, work, &
 134     2*(l/2+1)*m, ier1 )
 135
 136   if ( ier1 /= 0 ) then
 137      ier = 20
 138      call xerfft ( 'rfft2b', -5 )
 139      return
 140   end if
 141 !
 142 !  Reshuffle.
 143 !
 144   do j = 1, m
 145     do i = 2, l
 146       r(i,j) = r(i+1,j)
 147     end do
 148   end do
 149 !
 150 !  Transform first dimension of array.
 151 !
 152   call rfftmb ( m, ldim, l, 1, r, m*ldim, wsave(1), &
 153     l+int(log( real ( l, kind = 4 )))+4, work, 2*(l/2+1)*m, ier1 )
 154
 155   if ( ier1 /= 0 ) then
 156     ier = 20
 157     call xerfft ( 'rfft2f', -5 )
 158     return
 159   end if
 160
 161   return
 162 end