external/fftpack/fftpack5/sintb1.F

   1 subroutine sintb1 ( n, inc, x, wsave, xh, work, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 !! SINTB1 is an FFTPACK5 auxiliary routine.
   6 !
   7 !
   8 !    Copyright (C) 1995-2004, Scientific Computing Division,
   9 !    University Corporation for Atmospheric Research
  10 !
  11 !  Modified:
  12 !
  13 !    27 March 2009
  14 !
  15 !  Author:
  16 !
  17 !    Paul Swarztrauber
  18 !    Richard Valent
  19 !
  20 !  Reference:
  21 !
  22 !    Paul Swarztrauber,
  23 !    Vectorizing the Fast Fourier Transforms,
  24 !    in Parallel Computations,
  25 !    edited by G. Rodrigue,
  26 !    Academic Press, 1982.
  27 !
  28 !    Paul Swarztrauber,
  29 !    Fast Fourier Transform Algorithms for Vector Computers,
  30 !    Parallel Computing, pages 45-63, 1984.
  31 !
  32 !  Parameters:
  33 !
  34   implicit none
  35
  36   integer ( kind = 4 ) inc
  37
  38   real ( kind = 8 ) dsum
  39   real ( kind = 4 ) fnp1s4
  40   integer ( kind = 4 ) i
  41   integer ( kind = 4 ) ier
  42   integer ( kind = 4 ) ier1
  43   integer ( kind = 4 ) k
  44   integer ( kind = 4 ) kc
  45   integer ( kind = 4 ) lnsv
  46   integer ( kind = 4 ) lnwk
  47   integer ( kind = 4 ) lnxh
  48   integer ( kind = 4 ) modn
  49   integer ( kind = 4 ) n
  50   integer ( kind = 4 ) np1
  51   integer ( kind = 4 ) ns2
  52   real ( kind = 4 ) srt3s2
  53   real ( kind = 4 ) t1
  54   real ( kind = 4 ) t2
  55   real ( kind = 4 ) work(*)
  56   real ( kind = 4 ) wsave(*)
  57   real ( kind = 4 ) x(inc,*)
  58   real ( kind = 4 ) xh(*)
  59   real ( kind = 4 ) xhold
  60
  61   ier = 0
  62
  63   if ( n < 2 ) then
  64     return
  65   end if
  66
  67   if ( n == 2 ) then
  68     srt3s2 = sqrt ( 3.0E+00 ) / 2.0E+00
  69     xhold = srt3s2 * ( x(1,1) + x(1,2) )
  70     x(1,2) = srt3s2 * ( x(1,1) - x(1,2) )
  71     x(1,1) = xhold
  72     return
  73   end if
  74
  75   np1 = n + 1
  76   ns2 = n / 2
  77
  78   do k = 1, ns2
  79     kc = np1 - k
  80     t1 = x(1,k) - x(1,kc)
  81     t2 = wsave(k) * ( x(1,k) + x(1,kc) )
  82     xh(k+1) = t1 + t2
  83     xh(kc+1) = t2 - t1
  84   end do
  85
  86   modn = mod ( n, 2 )
  87
  88   if ( modn /= 0 ) then
  89     xh(ns2+2) = 4.0E+00 * x(1,ns2+1)
  90   end if
  91
  92   xh(1) = 0.0E+00
  93   lnxh = np1
  94   lnsv = np1 + int ( log ( real ( np1, kind = 4 ) ) ) + 4
  95   lnwk = np1
  96
  97   call rfft1f ( np1, 1, xh, lnxh, wsave(ns2+1), lnsv, work, lnwk, ier1 )
  98
  99   if ( ier1 /= 0 ) then
 100     ier = 20
 101     call xerfft ( 'sintb1', -5 )
 102     return
 103   end if
 104
 105   if ( mod ( np1, 2 ) == 0 ) then
 106     xh(np1) = xh(np1) + xh(np1)
 107   end if
 108
 109   fnp1s4 = real ( np1, kind = 4 ) / 4.0E+00
 110   x(1,1) = fnp1s4 * xh(1)
 111   dsum = x(1,1)
 112
 113   do i = 3, n, 2
 114     x(1,i-1) = fnp1s4 * xh(i)
 115     dsum = dsum + fnp1s4 * xh(i-1)
 116     x(1,i) = dsum
 117   end do
 118
 119   if ( modn == 0 ) then
 120     x(1,n) = fnp1s4 * xh(n+1)
 121   end if
 122
 123   return
 124 end