external/fftpack/fftpack5/r1fgkb.F

   1 subroutine r1fgkb ( ido, ip, l1, idl1, cc, c1, c2, in1, ch, ch2, in2, wa )
   2
   3 !*****************************************************************************80
   4 !
   5 !! R1FGKB 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 ) idl1
  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
  43   real ( kind = 4 ) ai1
  44   real ( kind = 4 ) ai2
  45   real ( kind = 4 ) ar1
  46   real ( kind = 4 ) ar1h
  47   real ( kind = 4 ) ar2
  48   real ( kind = 4 ) ar2h
  49   real ( kind = 4 ) arg
  50   real ( kind = 4 ) c1(in1,ido,l1,ip)
  51   real ( kind = 4 ) c2(in1,idl1,ip)
  52   real ( kind = 4 ) cc(in1,ido,ip,l1)
  53   real ( kind = 4 ) ch(in2,ido,l1,ip)
  54   real ( kind = 4 ) ch2(in2,idl1,ip)
  55   real ( kind = 4 ) dc2
  56   real ( kind = 4 ) dcp
  57   real ( kind = 4 ) ds2
  58   real ( kind = 4 ) dsp
  59   integer ( kind = 4 ) i
  60   integer ( kind = 4 ) ic
  61   integer ( kind = 4 ) idij
  62   integer ( kind = 4 ) idp2
  63   integer ( kind = 4 ) ik
  64   integer ( kind = 4 ) ipp2
  65   integer ( kind = 4 ) ipph
  66   integer ( kind = 4 ) is
  67   integer ( kind = 4 ) j
  68   integer ( kind = 4 ) j2
  69   integer ( kind = 4 ) jc
  70   integer ( kind = 4 ) k
  71   integer ( kind = 4 ) l
  72   integer ( kind = 4 ) lc
  73   integer ( kind = 4 ) nbd
  74   real ( kind = 4 ) tpi
  75   real ( kind = 4 ) wa(ido)
  76
  77   tpi = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 )
  78   arg = tpi / real ( ip, kind = 4 )
  79   dcp = cos ( arg )
  80   dsp = sin ( arg )
  81   idp2 = ido + 2
  82   nbd = ( ido - 1 ) / 2
  83   ipp2 = ip + 2
  84   ipph = ( ip + 1 ) / 2
  85
  86   if ( ido < l1 ) then
  87     do i = 1, ido
  88       do k = 1, l1
  89         ch(1,i,k,1) = cc(1,i,1,k)
  90       end do
  91     end do
  92   else
  93     do k = 1, l1
  94       do i = 1, ido
  95         ch(1,i,k,1) = cc(1,i,1,k)
  96       end do
  97     end do
  98   end if
  99
 100   do j = 2, ipph
 101     jc = ipp2 - j
 102     j2 = j + j
 103     do k = 1, l1
 104       ch(1,1,k,j) = cc(1,ido,j2-2,k)+cc(1,ido,j2-2,k)
 105       ch(1,1,k,jc) = cc(1,1,j2-1,k)+cc(1,1,j2-1,k)
 106     end do
 107   end do
 108
 109   if ( ido == 1 ) then
 110
 111   else if ( nbd < l1 ) then
 112
 113     do j = 2, ipph
 114       jc = ipp2 - j
 115       do i = 3, ido, 2
 116         ic = idp2 - i
 117         do k = 1, l1
 118           ch(1,i-1,k,j) = cc(1,i-1,2*j-1,k)+cc(1,ic-1,2*j-2,k)
 119           ch(1,i-1,k,jc) = cc(1,i-1,2*j-1,k)-cc(1,ic-1,2*j-2,k)
 120           ch(1,i,k,j) = cc(1,i,2*j-1,k)-cc(1,ic,2*j-2,k)
 121           ch(1,i,k,jc) = cc(1,i,2*j-1,k)+cc(1,ic,2*j-2,k)
 122         end do
 123       end do
 124     end do
 125
 126   else
 127
 128     do j = 2, ipph
 129       jc = ipp2 - j
 130       do k = 1, l1
 131         do i = 3, ido, 2
 132           ic = idp2 - i
 133           ch(1,i-1,k,j) = cc(1,i-1,2*j-1,k)+cc(1,ic-1,2*j-2,k)
 134           ch(1,i-1,k,jc) = cc(1,i-1,2*j-1,k)-cc(1,ic-1,2*j-2,k)
 135           ch(1,i,k,j) = cc(1,i,2*j-1,k)-cc(1,ic,2*j-2,k)
 136           ch(1,i,k,jc) = cc(1,i,2*j-1,k)+cc(1,ic,2*j-2,k)
 137         end do
 138       end do
 139     end do
 140
 141   end if
 142
 143   ar1 = 1.0E+00
 144   ai1 = 0.0E+00
 145
 146   do l = 2, ipph
 147
 148     lc = ipp2 - l
 149     ar1h = dcp * ar1 - dsp * ai1
 150     ai1 = dcp * ai1 + dsp * ar1
 151     ar1 = ar1h
 152
 153     do ik = 1, idl1
 154       c2(1,ik,l) = ch2(1,ik,1)+ar1*ch2(1,ik,2)
 155       c2(1,ik,lc) = ai1*ch2(1,ik,ip)
 156     end do
 157
 158     dc2 = ar1
 159     ds2 = ai1
 160     ar2 = ar1
 161     ai2 = ai1
 162
 163     do j = 3, ipph
 164
 165       jc = ipp2 - j
 166       ar2h = dc2*ar2-ds2*ai2
 167       ai2 = dc2*ai2+ds2*ar2
 168       ar2 = ar2h
 169
 170       do ik = 1, idl1
 171         c2(1,ik,l) = c2(1,ik,l)+ar2*ch2(1,ik,j)
 172         c2(1,ik,lc) = c2(1,ik,lc)+ai2*ch2(1,ik,jc)
 173       end do
 174
 175     end do
 176
 177   end do
 178
 179   do j = 2, ipph
 180     do ik = 1, idl1
 181       ch2(1,ik,1) = ch2(1,ik,1)+ch2(1,ik,j)
 182     end do
 183   end do
 184
 185   do j = 2, ipph
 186     jc = ipp2 - j
 187     do k = 1, l1
 188       ch(1,1,k,j) = c1(1,1,k,j)-c1(1,1,k,jc)
 189       ch(1,1,k,jc) = c1(1,1,k,j)+c1(1,1,k,jc)
 190     end do
 191   end do
 192
 193   if ( ido == 1 ) then
 194
 195   else if ( nbd < l1 ) then
 196
 197     do j = 2, ipph
 198       jc = ipp2 - j
 199       do i = 3, ido, 2
 200         do k = 1, l1
 201           ch(1,i-1,k,j)  = c1(1,i-1,k,j) - c1(1,i,k,jc)
 202           ch(1,i-1,k,jc) = c1(1,i-1,k,j) + c1(1,i,k,jc)
 203           ch(1,i,k,j)    = c1(1,i,k,j)   + c1(1,i-1,k,jc)
 204           ch(1,i,k,jc)   = c1(1,i,k,j)   - c1(1,i-1,k,jc)
 205         end do
 206       end do
 207     end do
 208
 209   else
 210
 211     do j = 2, ipph
 212       jc = ipp2 - j
 213       do k = 1, l1
 214         do i = 3, ido, 2
 215           ch(1,i-1,k,j) = c1(1,i-1,k,j)-c1(1,i,k,jc)
 216           ch(1,i-1,k,jc) = c1(1,i-1,k,j)+c1(1,i,k,jc)
 217           ch(1,i,k,j) = c1(1,i,k,j)+c1(1,i-1,k,jc)
 218           ch(1,i,k,jc) = c1(1,i,k,j)-c1(1,i-1,k,jc)
 219         end do
 220       end do
 221     end do
 222
 223   end if
 224
 225   if ( ido == 1 ) then
 226     return
 227   end if
 228
 229   do ik = 1, idl1
 230     c2(1,ik,1) = ch2(1,ik,1)
 231   end do
 232
 233   do j = 2, ip
 234     do k = 1, l1
 235       c1(1,1,k,j) = ch(1,1,k,j)
 236     end do
 237   end do
 238
 239   if ( l1 < nbd ) then
 240
 241     is = -ido
 242     do j = 2, ip
 243        is = is + ido
 244        do k = 1, l1
 245          idij = is
 246          do i = 3, ido, 2
 247            idij = idij + 2
 248            c1(1,i-1,k,j) = wa(idij-1)*ch(1,i-1,k,j)-wa(idij)* ch(1,i,k,j)
 249            c1(1,i,k,j) = wa(idij-1)*ch(1,i,k,j)+wa(idij)* ch(1,i-1,k,j)
 250          end do
 251        end do
 252     end do
 253
 254   else
 255
 256     is = -ido
 257
 258     do j = 2, ip
 259       is = is + ido
 260       idij = is
 261       do i = 3, ido, 2
 262         idij = idij + 2
 263         do k = 1, l1
 264            c1(1,i-1,k,j) = wa(idij-1) * ch(1,i-1,k,j) - wa(idij) * ch(1,i,k,j)
 265            c1(1,i,k,j)   = wa(idij-1) * ch(1,i,k,j)   + wa(idij) * ch(1,i-1,k,j)
 266         end do
 267       end do
 268     end do
 269
 270   end if
 271
 272   return
 273 end