external/fftpack/fftpack5/cosqf1.F

   1 subroutine cosqf1 ( n, inc, x, wsave, work, ier )
   2
   3 !*****************************************************************************80
   4 !
   5 !! COSQF1 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   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 ) lenx
  44   integer ( kind = 4 ) lnsv
  45   integer ( kind = 4 ) lnwk
  46   integer ( kind = 4 ) modn
  47   integer ( kind = 4 ) n
  48   integer ( kind = 4 ) np2
  49   integer ( kind = 4 ) ns2
  50   real ( kind = 4 ) work(*)
  51   real ( kind = 4 ) wsave(*)
  52   real ( kind = 4 ) x(inc,*)
  53   real ( kind = 4 ) xim1
  54
  55   ier = 0
  56   ns2 = ( n + 1 ) / 2
  57   np2 = n + 2
  58
  59   do k = 2, ns2
  60     kc = np2 - k
  61     work(k)  = x(1,k) + x(1,kc)
  62     work(kc) = x(1,k) - x(1,kc)
  63   end do
  64
  65   modn = mod ( n, 2 )
  66
  67   if ( modn == 0 ) then
  68     work(ns2+1) = x(1,ns2+1) + x(1,ns2+1)
  69   end if
  70
  71   do k = 2, ns2
  72     kc = np2 - k
  73     x(1,k)  = wsave(k-1) * work(kc) + wsave(kc-1) * work(k)
  74     x(1,kc) = wsave(k-1) * work(k)  - wsave(kc-1) * work(kc)
  75   end do
  76
  77   if ( modn == 0 ) then
  78     x(1,ns2+1) = wsave(ns2) * work(ns2+1)
  79   end if
  80
  81   lenx = inc * ( n - 1 ) + 1
  82   lnsv = n + int ( log ( real ( n, kind = 4 ) ) ) + 4
  83   lnwk = n
  84
  85   call rfft1f ( n, inc, x, lenx, wsave(n+1), lnsv, work, lnwk, ier1 )
  86
  87   if ( ier1 /= 0 ) then
  88     ier = 20
  89     call xerfft ( 'cosqf1', -5 )
  90     return
  91   end if
  92
  93   do i = 3, n, 2
  94     xim1   = 0.5E+00 * ( x(1,i-1) + x(1,i) )
  95     x(1,i) = 0.5E+00 * ( x(1,i-1) - x(1,i) )
  96     x(1,i-1) = xim1
  97   end do
  98
  99   return
 100 end