external/fftpack/fftpack5/dsintb1.F

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