external/fftpack/fftpack5/rfftb1.F

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