external/fftpack/fftpack5/dcosqf1.F

   1 subroutine dcosqf1 ( n, inc, x, wsave, work, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 !! DCOSQF1 is an FFTPACK5 auxiliary routine.
   6 !
   7 !
   8 !
   9 !  Modified:
  10 !
  11 !    17 November 2007
  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   integer ( kind = 4 ) i
  37   integer ( kind = 4 ) ier
  38   integer ( kind = 4 ) ier1
  39   integer ( kind = 4 ) k
  40   integer ( kind = 4 ) kc
  41   integer ( kind = 4 ) lenx
  42   integer ( kind = 4 ) lnsv
  43   integer ( kind = 4 ) lnwk
  44   integer ( kind = 4 ) modn
  45   integer ( kind = 4 ) n
  46   integer ( kind = 4 ) np2
  47   integer ( kind = 4 ) ns2
  48   real ( kind = 8 ) work(*)
  49   real ( kind = 8 ) wsave(*)
  50   real ( kind = 8 ) x(inc,*)
  51   real ( kind = 8 ) xim1
  52
  53   ier = 0
  54   ns2 = ( n + 1 ) / 2
  55   np2 = n + 2
  56
  57   do k = 2, ns2
  58     kc = np2 - k
  59     work(k)  = x(1,k) + x(1,kc)
  60     work(kc) = x(1,k) - x(1,kc)
  61   end do
  62
  63   modn = mod ( n, 2 )
  64
  65   if ( modn == 0 ) then
  66     work(ns2+1) = x(1,ns2+1) + x(1,ns2+1)
  67   end if
  68
  69   do k = 2, ns2
  70     kc = np2 - k
  71     x(1,k)  = wsave(k-1) * work(kc) + wsave(kc-1) * work(k)
  72     x(1,kc) = wsave(k-1) * work(k)  - wsave(kc-1) * work(kc)
  73   end do
  74
  75   if ( modn == 0 ) then
  76     x(1,ns2+1) = wsave(ns2) * work(ns2+1)
  77   end if
  78
  79   lenx = inc * ( n - 1 ) + 1
  80   lnsv = n + int ( log ( real ( n, kind = 8 ) ) ) + 4
  81   lnwk = n
  82
  83   call dfft1f ( n, inc, x, lenx, wsave(n+1), lnsv, work, lnwk, ier1 )
  84
  85   if ( ier1 /= 0 ) then
  86     ier = 20
  87     call xerfft ( 'dcosqf1', -5 )
  88     return
  89   end if
  90
  91   do i = 3, n, 2
  92     xim1   = 0.5D+00 * ( x(1,i-1) + x(1,i) )
  93     x(1,i) = 0.5D+00 * ( x(1,i-1) - x(1,i) )
  94     x(1,i-1) = xim1
  95   end do
  96
  97   return
  98 end