external/fftpack/fftpack5/zmfgkb.F

   1 subroutine zmfgkb ( lot, ido, ip, l1, lid, na, cc, cc1, im1, in1, &
   2   ch, ch1, im2, in2, wa )
   3
   4 !*****************************************************************************80
   5 !
   6 !! ZMFGKB is an FFTPACK5 auxiliary routine.
   7 !
   8 !
   9 !
  10 !  Modified:
  11 !
  12 !    26 Ausust 2009
  13 !
  14 !  Author:
  15 !
  16 !    Original complex single precision by Paul Swarztrauber, Richard Valent.
  17 !    Complex double precision version by John Burkardt.
  18 !
  19 !  Reference:
  20 !
  21 !    Paul Swarztrauber,
  22 !    Vectorizing the Fast Fourier Transforms,
  23 !    in Parallel Computations,
  24 !    edited by G. Rodrigue,
  25 !    Academic Press, 1982.
  26 !
  27 !    Paul Swarztrauber,
  28 !    Fast Fourier Transform Algorithms for Vector Computers,
  29 !    Parallel Computing, pages 45-63, 1984.
  30 !
  31 !  Parameters:
  32 !
  33   implicit none
  34
  35   integer ( kind = 4 ) ido
  36   integer ( kind = 4 ) in1
  37   integer ( kind = 4 ) in2
  38   integer ( kind = 4 ) ip
  39   integer ( kind = 4 ) l1
  40   integer ( kind = 4 ) lid
  41
  42   real ( kind = 8 ) cc(2,in1,l1,ip,ido)
  43   real ( kind = 8 ) cc1(2,in1,lid,ip)
  44   real ( kind = 8 ) ch(2,in2,l1,ido,ip)
  45   real ( kind = 8 ) ch1(2,in2,lid,ip)
  46   real ( kind = 8 ) chold1
  47   real ( kind = 8 ) chold2
  48   integer ( kind = 4 ) i
  49   integer ( kind = 4 ) idlj
  50   integer ( kind = 4 ) im1
  51   integer ( kind = 4 ) im2
  52   integer ( kind = 4 ) ipp2
  53   integer ( kind = 4 ) ipph
  54   integer ( kind = 4 ) j
  55   integer ( kind = 4 ) jc
  56   integer ( kind = 4 ) k
  57   integer ( kind = 4 ) ki
  58   integer ( kind = 4 ) l
  59   integer ( kind = 4 ) lc
  60   integer ( kind = 4 ) lot
  61   integer ( kind = 4 ) m1
  62   integer ( kind = 4 ) m1d
  63   integer ( kind = 4 ) m2
  64   integer ( kind = 4 ) m2s
  65   integer ( kind = 4 ) na
  66   real ( kind = 8 ) wa(ido,ip-1,2)
  67   real ( kind = 8 ) wai
  68   real ( kind = 8 ) war
  69
  70   m1d = ( lot - 1 ) * im1 + 1
  71   m2s = 1 - im2
  72   ipp2 = ip + 2
  73   ipph = ( ip + 1 ) / 2
  74
  75   do ki = 1, lid
  76     m2 = m2s
  77     do m1 = 1, m1d, im1
  78       m2 = m2 + im2
  79       ch1(1,m2,ki,1) = cc1(1,m1,ki,1)
  80       ch1(2,m2,ki,1) = cc1(2,m1,ki,1)
  81     end do
  82   end do
  83
  84   do j = 2, ipph
  85     jc = ipp2 - j
  86     do ki = 1, lid
  87       m2 = m2s
  88       do m1 = 1, m1d, im1
  89         m2 = m2 + im2
  90         ch1(1,m2,ki,j) =  cc1(1,m1,ki,j) + cc1(1,m1,ki,jc)
  91         ch1(1,m2,ki,jc) = cc1(1,m1,ki,j) - cc1(1,m1,ki,jc)
  92         ch1(2,m2,ki,j) =  cc1(2,m1,ki,j) + cc1(2,m1,ki,jc)
  93         ch1(2,m2,ki,jc) = cc1(2,m1,ki,j) - cc1(2,m1,ki,jc)
  94       end do
  95     end do
  96   end do
  97
  98   do j = 2, ipph
  99     do ki = 1, lid
 100       m2 = m2s
 101       do m1 = 1, m1d, im1
 102         m2 = m2 + im2
 103         cc1(1,m1,ki,1) = cc1(1,m1,ki,1) + ch1(1,m2,ki,j)
 104         cc1(2,m1,ki,1) = cc1(2,m1,ki,1) + ch1(2,m2,ki,j)
 105       end do
 106     end do
 107   end do
 108
 109   do l = 2, ipph
 110
 111     lc = ipp2 - l
 112     do ki = 1, lid
 113       m2 = m2s
 114       do m1 = 1, m1d, im1
 115         m2 = m2 + im2
 116
 117         cc1(1,m1,ki,l)  = ch1(1,m2,ki,1) + wa(1,l-1,1) * ch1(1,m2,ki,2)
 118         cc1(1,m1,ki,lc) =                  wa(1,l-1,2) * ch1(1,m2,ki,ip)
 119         cc1(2,m1,ki,l)  = ch1(2,m2,ki,1) + wa(1,l-1,1) * ch1(2,m2,ki,2)
 120         cc1(2,m1,ki,lc) =                  wa(1,l-1,2) * ch1(2,m2,ki,ip)
 121
 122       end do
 123     end do
 124
 125     do j = 3, ipph
 126       jc = ipp2 - j
 127       idlj = mod ( ( l - 1 ) * ( j - 1 ), ip )
 128       war = wa(1,idlj,1)
 129       wai = wa(1,idlj,2)
 130       do ki = 1, lid
 131         m2 = m2s
 132         do m1 = 1, m1d, im1
 133           m2 = m2 + im2
 134           cc1(1,m1,ki,l)  = cc1(1,m1,ki,l)  + war * ch1(1,m2,ki,j)
 135           cc1(1,m1,ki,lc) = cc1(1,m1,ki,lc) + wai * ch1(1,m2,ki,jc)
 136           cc1(2,m1,ki,l)  = cc1(2,m1,ki,l)  + war * ch1(2,m2,ki,j)
 137           cc1(2,m1,ki,lc) = cc1(2,m1,ki,lc) + wai * ch1(2,m2,ki,jc)
 138         end do
 139       end do
 140     end do
 141
 142   end do
 143
 144   if( 1 < ido .or. na == 1 ) then
 145
 146     do ki = 1, lid
 147       m2 = m2s
 148       do m1 = 1, m1d, im1
 149         m2 = m2 + im2
 150         ch1(1,m2,ki,1) = cc1(1,m1,ki,1)
 151         ch1(2,m2,ki,1) = cc1(2,m1,ki,1)
 152       end do
 153     end do
 154
 155     do j = 2, ipph
 156       jc = ipp2 - j
 157       do ki = 1, lid
 158         m2 = m2s
 159         do m1 = 1, m1d, im1
 160           m2 = m2 + im2
 161           ch1(1,m2,ki,j)  = cc1(1,m1,ki,j) - cc1(2,m1,ki,jc)
 162           ch1(1,m2,ki,jc) = cc1(1,m1,ki,j) + cc1(2,m1,ki,jc)
 163           ch1(2,m2,ki,jc) = cc1(2,m1,ki,j) - cc1(1,m1,ki,jc)
 164           ch1(2,m2,ki,j)  = cc1(2,m1,ki,j) + cc1(1,m1,ki,jc)
 165         end do
 166       end do
 167     end do
 168
 169     if ( ido == 1 ) then
 170       return
 171     end if
 172
 173     do i = 1, ido
 174       do k = 1, l1
 175         m2 = m2s
 176         do m1 = 1, m1d, im1
 177           m2 = m2 + im2
 178           cc(1,m1,k,1,i) = ch(1,m2,k,i,1)
 179           cc(2,m1,k,1,i) = ch(2,m2,k,i,1)
 180         end do
 181       end do
 182     end do
 183
 184     do j = 2, ip
 185       do k = 1, l1
 186         m2 = m2s
 187         do m1 = 1, m1d, im1
 188           m2 = m2 + im2
 189           cc(1,m1,k,j,1) = ch(1,m2,k,1,j)
 190           cc(2,m1,k,j,1) = ch(2,m2,k,1,j)
 191         end do
 192       end do
 193     end do
 194
 195     do j = 2, ip
 196       do i = 2, ido
 197         do k = 1, l1
 198           m2 = m2s
 199           do m1 = 1, m1d, im1
 200             m2 = m2 + im2
 201             cc(1,m1,k,j,i) = wa(i,j-1,1) * ch(1,m2,k,i,j) &
 202                            - wa(i,j-1,2) * ch(2,m2,k,i,j)
 203             cc(2,m1,k,j,i) = wa(i,j-1,1) * ch(2,m2,k,i,j) &
 204                            + wa(i,j-1,2) * ch(1,m2,k,i,j)
 205           end do
 206         end do
 207       end do
 208     end do
 209
 210   else
 211
 212     do j = 2, ipph
 213       jc = ipp2 - j
 214       do ki = 1, lid
 215         do m1 = 1, m1d, im1
 216
 217           chold1         = cc1(1,m1,ki,j) - cc1(2,m1,ki,jc)
 218           chold2         = cc1(1,m1,ki,j) + cc1(2,m1,ki,jc)
 219           cc1(1,m1,ki,j) = chold1
 220
 221           cc1(2,m1,ki,jc) = cc1(2,m1,ki,j) - cc1(1,m1,ki,jc)
 222           cc1(2,m1,ki,j)  = cc1(2,m1,ki,j) + cc1(1,m1,ki,jc)
 223           cc1(1,m1,ki,jc) = chold2
 224
 225         end do
 226       end do
 227     end do
 228
 229   end if
 230
 231   return
 232 end