external/fftpack/fftpack5/mrftb1.F

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