external/fftpack/fftpack5/dfftb1.F

   1 subroutine dfftb1 ( n, in, c, ch, wa, fac )
   2
   3 !*****************************************************************************80
   4 !
   5 !! DFFTB1 is an FFTPACK5 auxiliary routine.
   6 !
   7 !
   8 !
   9 !  Modified:
  10 !
  11 !    27 March 2009
  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   real ( kind = 8 ) half
  41   real ( kind = 8 ) halfm
  42   integer ( kind = 4 ) idl1
  43   integer ( kind = 4 ) ido
  44   integer ( kind = 4 ) ip
  45   integer ( kind = 4 ) iw
  46   integer ( kind = 4 ) ix2
  47   integer ( kind = 4 ) ix3
  48   integer ( kind = 4 ) ix4
  49   integer ( kind = 4 ) j
  50   integer ( kind = 4 ) k1
  51   integer ( kind = 4 ) l1
  52   integer ( kind = 4 ) l2
  53   integer ( kind = 4 ) modn
  54   integer ( kind = 4 ) na
  55   integer ( kind = 4 ) nf
  56   integer ( kind = 4 ) nl
  57   real ( kind = 8 ) wa(n)
  58
  59   nf = int ( fac(2) )
  60   na = 0
  61
  62   do k1 = 1, nf
  63
  64     ip = int ( fac(k1+2) )
  65     na = 1 - na
  66
  67     if ( 5 < ip ) then
  68       if ( k1 /= nf ) then
  69         na = 1 - na
  70       end if
  71     end if
  72
  73   end do
  74
  75   half = 0.5D+00
  76   halfm = -0.5D+00
  77   modn = mod ( n, 2 )
  78   nl = n - 2
  79   if ( modn /= 0 ) then
  80     nl = n - 1
  81   end if
  82
  83   if ( na == 0 ) then
  84
  85     do j = 2, nl, 2
  86       c(1,j) = half * c(1,j)
  87       c(1,j+1) = halfm * c(1,j+1)
  88     end do
  89
  90   else
  91
  92     ch(1) = c(1,1)
  93     ch(n) = c(1,n)
  94
  95     do j = 2, nl, 2
  96       ch(j) = half * c(1,j)
  97       ch(j+1) = halfm * c(1,j+1)
  98     end do
  99
 100   end if
 101
 102   l1 = 1
 103   iw = 1
 104
 105   do k1 = 1, nf
 106
 107     ip = int ( fac(k1+2) )
 108     l2 = ip * l1
 109     ido = n / l2
 110     idl1 = ido * l1
 111
 112     if ( ip == 4 ) then
 113
 114       ix2 = iw + ido
 115       ix3 = ix2 + ido
 116
 117       if ( na == 0 ) then
 118         call d1f4kb ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2), wa(ix3) )
 119       else
 120         call d1f4kb ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2), wa(ix3) )
 121       end if
 122
 123       na = 1 - na
 124
 125     else if ( ip == 2 ) then
 126
 127       if ( na == 0 ) then
 128         call d1f2kb ( ido, l1, c, in, ch, 1, wa(iw) )
 129       else
 130         call d1f2kb ( ido, l1, ch, 1, c, in, wa(iw) )
 131       end if
 132
 133       na = 1 - na
 134
 135     else if ( ip == 3 ) then
 136
 137       ix2 = iw + ido
 138
 139       if ( na == 0 ) then
 140         call d1f3kb ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2) )
 141       else
 142         call d1f3kb ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2) )
 143       end if
 144
 145       na = 1 - na
 146
 147     else if ( ip == 5 ) then
 148
 149       ix2 = iw + ido
 150       ix3 = ix2 + ido
 151       ix4 = ix3 + ido
 152
 153       if ( na == 0 ) then
 154         call d1f5kb ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2), wa(ix3), wa(ix4) )
 155       else
 156         call d1f5kb ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2), wa(ix3), wa(ix4) )
 157       end if
 158
 159       na = 1 - na
 160
 161     else
 162
 163       if ( na == 0 ) then
 164         call d1fgkb ( ido, ip, l1, idl1, c, c, c, in, ch, ch, 1, wa(iw) )
 165       else
 166         call d1fgkb ( ido, ip, l1, idl1, ch, ch, ch, 1, c, c, in, wa(iw) )
 167       end if
 168
 169       if ( ido == 1 ) then
 170         na = 1 - na
 171       end if
 172
 173     end if
 174
 175     l1 = l2
 176     iw = iw + ( ip - 1 ) * ido
 177
 178   end do
 179
 180   return
 181 end