external/fftpack/fftpack5/zmfgkf.F

   1 subroutine zmfgkf ( lot, ido, ip, l1, lid, na, cc, cc1, im1, in1, &
   2   ch, ch1, im2, in2, wa )
   3
   4 !*****************************************************************************80
   5 !
   6 !! ZMFGKF 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 ) sn
  67   real ( kind = 8 ) wa(ido,ip-1,2)
  68   real ( kind = 8 ) wai
  69   real ( kind = 8 ) war
  70
  71   m1d = ( lot - 1 ) * im1 + 1
  72   m2s = 1 - im2
  73   ipp2 = ip + 2
  74   ipph = ( ip + 1 ) / 2
  75
  76   do ki = 1, lid
  77     m2 = m2s
  78     do m1 = 1, m1d, im1
  79       m2 = m2 + im2
  80       ch1(1,m2,ki,1) = cc1(1,m1,ki,1)
  81       ch1(2,m2,ki,1) = cc1(2,m1,ki,1)
  82     end do
  83   end do
  84
  85   do j = 2, ipph
  86     jc = ipp2 - j
  87     do ki = 1, lid
  88       m2 = m2s
  89       do m1 = 1, m1d, im1
  90         m2 = m2 + im2
  91         ch1(1,m2,ki,j) =  cc1(1,m1,ki,j) + cc1(1,m1,ki,jc)
  92         ch1(1,m2,ki,jc) = cc1(1,m1,ki,j) - cc1(1,m1,ki,jc)
  93         ch1(2,m2,ki,j) =  cc1(2,m1,ki,j) + cc1(2,m1,ki,jc)
  94         ch1(2,m2,ki,jc) = cc1(2,m1,ki,j) - cc1(2,m1,ki,jc)
  95       end do
  96     end do
  97   end do
  98
  99   do j = 2, ipph
 100     do ki = 1, lid
 101       m2 = m2s
 102       do m1 = 1, m1d, im1
 103         m2 = m2 + im2
 104         cc1(1,m1,ki,1) = cc1(1,m1,ki,1) + ch1(1,m2,ki,j)
 105         cc1(2,m1,ki,1) = cc1(2,m1,ki,1) + ch1(2,m2,ki,j)
 106       end do
 107     end do
 108   end do
 109
 110   do l = 2, ipph
 111
 112     lc = ipp2 - l
 113
 114     do ki = 1, lid
 115       m2 = m2s
 116       do m1 = 1, m1d, im1
 117         m2 = m2 + im2
 118         cc1(1,m1,ki,l)  = ch1(1,m2,ki,1) + wa(1,l-1,1) * ch1(1,m2,ki,2)
 119         cc1(1,m1,ki,lc) =                - wa(1,l-1,2) * ch1(1,m2,ki,ip)
 120         cc1(2,m1,ki,l)  = ch1(2,m2,ki,1) + wa(1,l-1,1) * ch1(2,m2,ki,2)
 121         cc1(2,m1,ki,lc) =                - wa(1,l-1,2) * ch1(2,m2,ki,ip)
 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 ) 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(2,m2,ki,j)  = cc1(2,m1,ki,j) + cc1(1,m1,ki,jc)
 163           ch1(1,m2,ki,jc) = cc1(1,m1,ki,j) + cc1(2,m1,ki,jc)
 164           ch1(2,m2,ki,jc) = cc1(2,m1,ki,j) - cc1(1,m1,ki,jc)
 165         end do
 166       end do
 167     end do
 168
 169     do i = 1, ido
 170       do k = 1, l1
 171         m2 = m2s
 172         do m1 = 1, m1d, im1
 173           m2 = m2 + im2
 174           cc(1,m1,k,1,i) = ch(1,m2,k,i,1)
 175           cc(2,m1,k,1,i) = ch(2,m2,k,i,1)
 176         end do
 177       end do
 178     end do
 179
 180     do j = 2, ip
 181       do k = 1, l1
 182         m2 = m2s
 183         do m1 = 1, m1d, im1
 184           m2 = m2 + im2
 185           cc(1,m1,k,j,1) = ch(1,m2,k,1,j)
 186           cc(2,m1,k,j,1) = ch(2,m2,k,1,j)
 187         end do
 188       end do
 189     end do
 190
 191     do j = 2, ip
 192       do i = 2, ido
 193         do k = 1, l1
 194           m2 = m2s
 195           do m1 = 1, m1d, im1
 196             m2 = m2 + im2
 197             cc(1,m1,k,j,i) = wa(i,j-1,1) * ch(1,m2,k,i,j) &
 198                            + wa(i,j-1,2) * ch(2,m2,k,i,j)
 199             cc(2,m1,k,j,i) = wa(i,j-1,1) * ch(2,m2,k,i,j) &
 200                            - wa(i,j-1,2) * ch(1,m2,k,i,j)
 201           end do
 202         end do
 203       end do
 204     end do
 205
 206   else if ( na == 1 ) then
 207
 208     sn = 1.0D+00 / real ( ip * l1, kind = 8 )
 209
 210     do ki = 1, lid
 211       m2 = m2s
 212       do m1 = 1, m1d, im1
 213         m2 = m2 + im2
 214         ch1(1,m2,ki,1) = sn * cc1(1,m1,ki,1)
 215         ch1(2,m2,ki,1) = sn * cc1(2,m1,ki,1)
 216       end do
 217     end do
 218
 219     do j = 2, ipph
 220       jc = ipp2 - j
 221       do ki = 1, lid
 222         m2 = m2s
 223         do m1 = 1, m1d, im1
 224           m2 = m2 + im2
 225           ch1(1,m2,ki,j)  = sn * ( cc1(1,m1,ki,j) - cc1(2,m1,ki,jc) )
 226           ch1(2,m2,ki,j)  = sn * ( cc1(2,m1,ki,j) + cc1(1,m1,ki,jc) )
 227           ch1(1,m2,ki,jc) = sn * ( cc1(1,m1,ki,j) + cc1(2,m1,ki,jc) )
 228           ch1(2,m2,ki,jc) = sn * ( cc1(2,m1,ki,j) - cc1(1,m1,ki,jc) )
 229         end do
 230       end do
 231     end do
 232
 233   else
 234
 235     sn = 1.0D+00 / real ( ip * l1, kind = 8 )
 236
 237     do ki = 1, lid
 238       m2 = m2s
 239       do m1 = 1, m1d, im1
 240         m2 = m2 + im2
 241         cc1(1,m1,ki,1) = sn * cc1(1,m1,ki,1)
 242         cc1(2,m1,ki,1) = sn * cc1(2,m1,ki,1)
 243       end do
 244     end do
 245
 246     do j = 2, ipph
 247       jc = ipp2 - j
 248       do ki = 1, lid
 249         do m1 = 1, m1d, im1
 250           chold1 = sn * ( cc1(1,m1,ki,j) - cc1(2,m1,ki,jc) )
 251           chold2 = sn * ( cc1(1,m1,ki,j) + cc1(2,m1,ki,jc) )
 252           cc1(1,m1,ki,j) = chold1
 253           cc1(2,m1,ki,jc) = sn * ( cc1(2,m1,ki,j) - cc1(1,m1,ki,jc) )
 254           cc1(2,m1,ki,j)  = sn * ( cc1(2,m1,ki,j) + cc1(1,m1,ki,jc) )
 255           cc1(1,m1,ki,jc) = chold2
 256         end do
 257       end do
 258     end do
 259
 260   end if
 261
 262   return
 263 end