external/fftpack/fftpack5/dcostf1.F

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