external/fftpack/fftpack5/r1fgkf.F

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