external/fftpack/fftpack5/cmfgkb.F

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