external/fftpack/fftpack5/dcostb1.F

   1 subroutine dcostb1 ( n, inc, x, wsave, work, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 !! DCOSTB1 is an FFTPACK5 auxiliary routine.
   6 !
   7 !
   8 !
   9 !  Modified:
  10 !
  11 !    07 February 2006
  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 ) fnm1s2
  38   real ( kind = 8 ) fnm1s4
  39   integer ( kind = 4 ) i
  40   integer ( kind = 4 ) ier
  41   integer ( kind = 4 ) ier1
  42   integer ( kind = 4 ) k
  43   integer ( kind = 4 ) kc
  44   integer ( kind = 4 ) lenx
  45   integer ( kind = 4 ) lnsv
  46   integer ( kind = 4 ) lnwk
  47   integer ( kind = 4 ) modn
  48   integer ( kind = 4 ) n
  49   integer ( kind = 4 ) nm1
  50   integer ( kind = 4 ) np1
  51   integer ( kind = 4 ) ns2
  52   real ( kind = 8 ) t1
  53   real ( kind = 8 ) t2
  54   real ( kind = 8 ) work(*)
  55   real ( kind = 8 ) wsave(*)
  56   real ( kind = 8 ) x(inc,*)
  57   real ( kind = 8 ) x1h
  58   real ( kind = 8 ) x1p3
  59   real ( kind = 8 ) x2
  60   real ( kind = 8 ) xi
  61
  62   ier = 0
  63   nm1 = n - 1
  64   np1 = n + 1
  65   ns2 = n / 2
  66
  67   if ( n < 2 ) then
  68     return
  69   end if
  70
  71   if ( n == 2 ) then
  72     x1h    = x(1,1) + x(1,2)
  73     x(1,2) = x(1,1) - x(1,2)
  74     x(1,1) = x1h
  75     return
  76   end if
  77
  78   if ( n == 3 ) then
  79     x1p3 = x(1,1) + x(1,3)
  80     x2 = x(1,2)
  81     x(1,2) = x(1,1) - x(1,3)
  82     x(1,1) = x1p3 + x2
  83     x(1,3) = x1p3 - x2
  84     return
  85   end if
  86
  87   x(1,1) = x(1,1) + x(1,1)
  88   x(1,n) = x(1,n) + x(1,n)
  89   dsum = x(1,1) - x(1,n)
  90   x(1,1) = x(1,1) + x(1,n)
  91
  92   do k = 2, ns2
  93     kc = np1 - k
  94     t1 = x(1,k) + x(1,kc)
  95     t2 = x(1,k) - x(1,kc)
  96     dsum = dsum + wsave(kc) * t2
  97     t2 = wsave(k) * t2
  98     x(1,k) = t1 - t2
  99     x(1,kc) = t1 + t2
 100   end do
 101
 102   modn = mod ( n, 2 )
 103
 104   if ( modn /= 0 ) then
 105     x(1,ns2+1) = x(1,ns2+1) + x(1,ns2+1)
 106   end if
 107
 108   lenx = inc * ( nm1 - 1 )  + 1
 109   lnsv = nm1 + int ( log ( real ( nm1, kind = 8 ) ) ) + 4
 110   lnwk = nm1
 111
 112   call dfft1f ( nm1, inc, x, lenx, wsave(n+1), lnsv, work, lnwk, ier1 )
 113
 114   if ( ier1 /= 0 ) then
 115     ier = 20
 116     call xerfft ( 'DCOSTB1', -5 )
 117     return
 118   end if
 119
 120   fnm1s2 = real ( nm1, kind = 8 ) / 2.0D+00
 121   dsum = 0.5D+00 * dsum
 122   x(1,1) = fnm1s2 * x(1,1)
 123
 124   if ( mod ( nm1, 2 ) == 0 ) then
 125     x(1,nm1) = x(1,nm1) + x(1,nm1)
 126   end if
 127
 128   fnm1s4 = real ( nm1, kind = 8 ) / 4.0D+00
 129
 130   do i = 3, n, 2
 131     xi = fnm1s4 * x(1,i)
 132     x(1,i) = fnm1s4 * x(1,i-1)
 133     x(1,i-1) = dsum
 134     dsum = dsum + xi
 135   end do
 136
 137   if ( modn == 0 ) then
 138     x(1,n) = dsum
 139   end if
 140
 141   return
 142 end