external/fftpack/fftpack5/dfftf1.F

   1 subroutine dfftf1 ( n, in, c, ch, wa, fac )
   2
   3 !*****************************************************************************80
   4 !
   5 !! DFFTF1 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 ) in
  35   integer ( kind = 4 ) n
  36
  37   real ( kind = 8 ) c(in,*)
  38   real ( kind = 8 ) ch(*)
  39   real ( kind = 8 ) fac(15)
  40   integer ( kind = 4 ) idl1
  41   integer ( kind = 4 ) ido
  42   integer ( kind = 4 ) ip
  43   integer ( kind = 4 ) iw
  44   integer ( kind = 4 ) ix2
  45   integer ( kind = 4 ) ix3
  46   integer ( kind = 4 ) ix4
  47   integer ( kind = 4 ) j
  48   integer ( kind = 4 ) k1
  49   integer ( kind = 4 ) kh
  50   integer ( kind = 4 ) l1
  51   integer ( kind = 4 ) l2
  52   integer ( kind = 4 ) modn
  53   integer ( kind = 4 ) na
  54   integer ( kind = 4 ) nf
  55   integer ( kind = 4 ) nl
  56   real ( kind = 8 ) sn
  57   real ( kind = 8 ) tsn
  58   real ( kind = 8 ) tsnm
  59   real ( kind = 8 ) wa(n)
  60
  61   nf = int ( fac(2) )
  62   na = 1
  63   l2 = n
  64   iw = n
  65
  66   do k1 = 1, nf
  67
  68     kh = nf - k1
  69     ip = int ( fac(kh+3) )
  70     l1 = l2 / ip
  71     ido = n / l2
  72     idl1 = ido * l1
  73     iw = iw - ( ip - 1 ) * ido
  74     na = 1 - na
  75
  76     if ( ip == 4 ) then
  77
  78       ix2 = iw + ido
  79       ix3 = ix2 + ido
  80
  81       if ( na == 0 ) then
  82         call d1f4kf ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2), wa(ix3) )
  83       else
  84         call d1f4kf ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2), wa(ix3) )
  85       end if
  86
  87     else if ( ip == 2 ) then
  88
  89       if ( na == 0 ) then
  90         call d1f2kf ( ido, l1, c, in, ch, 1, wa(iw) )
  91       else
  92         call d1f2kf ( ido, l1, ch, 1, c, in, wa(iw) )
  93       end if
  94
  95     else if ( ip == 3 ) then
  96
  97       ix2 = iw + ido
  98
  99       if ( na == 0 ) then
 100         call d1f3kf ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2) )
 101       else
 102         call d1f3kf ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2) )
 103       end if
 104
 105     else if ( ip == 5 ) then
 106
 107       ix2 = iw + ido
 108       ix3 = ix2 + ido
 109       ix4 = ix3 + ido
 110
 111       if ( na == 0 ) then
 112         call d1f5kf ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2), wa(ix3), wa(ix4) )
 113       else
 114         call d1f5kf ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2), wa(ix3), wa(ix4) )
 115       end if
 116
 117     else
 118
 119       if ( ido == 1 ) then
 120         na = 1 - na
 121       end if
 122
 123       if ( na == 0 ) then
 124         call d1fgkf ( ido, ip, l1, idl1, c, c, c, in, ch, ch, 1, wa(iw) )
 125         na = 1
 126       else
 127         call d1fgkf ( ido, ip, l1, idl1, ch, ch, ch, 1, c, c, in, wa(iw) )
 128         na = 0
 129       end if
 130
 131     end if
 132
 133     l2 = l1
 134
 135   end do
 136
 137   sn = 1.0D+00 / real ( n, kind = 8 )
 138   tsn = 2.0D+00 / real ( n, kind = 8 )
 139   tsnm = - tsn
 140   modn = mod ( n, 2 )
 141   nl = n - 2
 142
 143   if ( modn /= 0 ) then
 144     nl = n - 1
 145   end if
 146
 147   if ( na == 0 ) then
 148
 149     c(1,1) = sn * ch(1)
 150     do j = 2, nl, 2
 151       c(1,j) = tsn * ch(j)
 152       c(1,j+1) = tsnm * ch(j+1)
 153     end do
 154
 155     if ( modn == 0 ) then
 156       c(1,n) = sn * ch(n)
 157     end if
 158
 159   else
 160
 161     c(1,1) = sn * c(1,1)
 162
 163     do j = 2, nl, 2
 164       c(1,j) = tsn * c(1,j)
 165       c(1,j+1) = tsnm * c(1,j+1)
 166     end do
 167
 168     if ( modn == 0 ) then
 169       c(1,n) = sn * c(1,n)
 170     end if
 171
 172   end if
 173
 174   return
 175 end