external/fftpack/fftpack5/d1fgkf.F

   1 subroutine d1fgkf ( ido, ip, l1, idl1, cc, c1, c2, in1, ch, ch2, in2, wa )
   2
   3 !*****************************************************************************80
   4 !
   5 !! D1FGKF is an FFTPACK5 auxiliary routine.
   6 !
   7 !
   8 !
   9 !  Modified:
  10 !
  11 !    07 February 2006
  12 !
  13 !  Author:
  14 !
  15 !    Original real single precision by Paul Swarztrauber, Richard Valent.
  16 !    Real double precision version by John Burkardt.
  17 !
  18 !  Reference:
  19 !
  20 !    Paul Swarztrauber,
  21 !    Vectorizing the Fast Fourier Transforms,
  22 !    in Parallel Computations,
  23 !    edited by G. Rodrigue,
  24 !    Academic Press, 1982.
  25 !
  26 !    Paul Swarztrauber,
  27 !    Fast Fourier Transform Algorithms for Vector Computers,
  28 !    Parallel Computing, pages 45-63, 1984.
  29 !
  30 !  Parameters:
  31 !
  32   implicit none
  33
  34   integer ( kind = 4 ) idl1
  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
  41   real ( kind = 8 ) ai1
  42   real ( kind = 8 ) ai2
  43   real ( kind = 8 ) ar1
  44   real ( kind = 8 ) ar1h
  45   real ( kind = 8 ) ar2
  46   real ( kind = 8 ) ar2h
  47   real ( kind = 8 ) arg
  48   real ( kind = 8 ) c1(in1,ido,l1,ip)
  49   real ( kind = 8 ) c2(in1,idl1,ip)
  50   real ( kind = 8 ) cc(in1,ido,ip,l1)
  51   real ( kind = 8 ) ch(in2,ido,l1,ip)
  52   real ( kind = 8 ) ch2(in2,idl1,ip)
  53   real ( kind = 8 ) dc2
  54   real ( kind = 8 ) dcp
  55   real ( kind = 8 ) ds2
  56   real ( kind = 8 ) dsp
  57   integer ( kind = 4 ) i
  58   integer ( kind = 4 ) ic
  59   integer ( kind = 4 ) idij
  60   integer ( kind = 4 ) idp2
  61   integer ( kind = 4 ) ik
  62   integer ( kind = 4 ) ipp2
  63   integer ( kind = 4 ) ipph
  64   integer ( kind = 4 ) is
  65   integer ( kind = 4 ) j
  66   integer ( kind = 4 ) j2
  67   integer ( kind = 4 ) jc
  68   integer ( kind = 4 ) k
  69   integer ( kind = 4 ) l
  70   integer ( kind = 4 ) lc
  71   integer ( kind = 4 ) nbd
  72   real ( kind = 8 ) tpi
  73   real ( kind = 8 ) wa(ido)
  74
  75   tpi = 2.0D+00 * 4.0D+00 * atan ( 1.0D+00 )
  76   arg = tpi / real ( ip, kind = 8 )
  77   dcp = cos ( arg )
  78   dsp = sin ( arg )
  79   ipph = ( ip + 1 ) / 2
  80   ipp2 = ip + 2
  81   idp2 = ido + 2
  82   nbd = ( ido - 1 ) / 2
  83
  84   if ( ido == 1 ) then
  85
  86     do ik = 1, idl1
  87       c2(1,ik,1) = ch2(1,ik,1)
  88     end do
  89
  90   else
  91
  92     do ik = 1, idl1
  93       ch2(1,ik,1) = c2(1,ik,1)
  94     end do
  95
  96     do j = 2, ip
  97       do k = 1, l1
  98         ch(1,1,k,j) = c1(1,1,k,j)
  99       end do
 100     end do
 101
 102     if ( l1 < nbd ) then
 103
 104       is = -ido
 105
 106       do j = 2, ip
 107         is = is + ido
 108         do k = 1, l1
 109           idij = is
 110           do i = 3, ido, 2
 111             idij = idij + 2
 112             ch(1,i-1,k,j) = wa(idij-1)*c1(1,i-1,k,j)+wa(idij) *c1(1,i,k,j)
 113             ch(1,i,k,j) = wa(idij-1)*c1(1,i,k,j)-wa(idij) *c1(1,i-1,k,j)
 114           end do
 115         end do
 116       end do
 117
 118     else
 119
 120       is = -ido
 121
 122       do j = 2, ip
 123         is = is + ido
 124         idij = is
 125         do i = 3, ido, 2
 126           idij = idij + 2
 127           do k = 1, l1
 128             ch(1,i-1,k,j) = wa(idij-1)*c1(1,i-1,k,j)+wa(idij) *c1(1,i,k,j)
 129             ch(1,i,k,j) = wa(idij-1)*c1(1,i,k,j)-wa(idij) *c1(1,i-1,k,j)
 130           end do
 131         end do
 132       end do
 133
 134     end if
 135
 136     if ( nbd < l1 ) then
 137
 138       do j = 2, ipph
 139         jc = ipp2 - j
 140         do i = 3, ido, 2
 141           do k = 1, l1
 142             c1(1,i-1,k,j) = ch(1,i-1,k,j)+ch(1,i-1,k,jc)
 143             c1(1,i-1,k,jc) = ch(1,i,k,j)-ch(1,i,k,jc)
 144             c1(1,i,k,j) = ch(1,i,k,j)+ch(1,i,k,jc)
 145             c1(1,i,k,jc) = ch(1,i-1,k,jc)-ch(1,i-1,k,j)
 146           end do
 147         end do
 148       end do
 149
 150     else
 151
 152       do j = 2, ipph
 153         jc = ipp2 - j
 154         do k = 1, l1
 155           do i = 3, ido, 2
 156             c1(1,i-1,k,j) = ch(1,i-1,k,j)+ch(1,i-1,k,jc)
 157             c1(1,i-1,k,jc) = ch(1,i,k,j)-ch(1,i,k,jc)
 158             c1(1,i,k,j) = ch(1,i,k,j)+ch(1,i,k,jc)
 159             c1(1,i,k,jc) = ch(1,i-1,k,jc)-ch(1,i-1,k,j)
 160           end do
 161         end do
 162       end do
 163
 164     end if
 165
 166   end if
 167
 168   do j = 2, ipph
 169     jc = ipp2 - j
 170     do k = 1, l1
 171       c1(1,1,k,j) = ch(1,1,k,j)+ch(1,1,k,jc)
 172       c1(1,1,k,jc) = ch(1,1,k,jc)-ch(1,1,k,j)
 173     end do
 174   end do
 175
 176   ar1 = 1.0D+00
 177   ai1 = 0.0D+00
 178
 179   do l = 2, ipph
 180
 181     lc = ipp2 - l
 182     ar1h = dcp * ar1 - dsp * ai1
 183     ai1 = dcp * ai1 + dsp * ar1
 184     ar1 = ar1h
 185
 186     do ik = 1, idl1
 187       ch2(1,ik,l) = c2(1,ik,1)+ar1*c2(1,ik,2)
 188       ch2(1,ik,lc) = ai1*c2(1,ik,ip)
 189     end do
 190
 191     dc2 = ar1
 192     ds2 = ai1
 193     ar2 = ar1
 194     ai2 = ai1
 195
 196     do j = 3, ipph
 197       jc = ipp2 - j
 198       ar2h = dc2 * ar2 - ds2 * ai2
 199       ai2 = dc2 * ai2 + ds2 * ar2
 200       ar2 = ar2h
 201       do ik = 1, idl1
 202         ch2(1,ik,l) = ch2(1,ik,l)+ar2*c2(1,ik,j)
 203         ch2(1,ik,lc) = ch2(1,ik,lc)+ai2*c2(1,ik,jc)
 204       end do
 205     end do
 206
 207   end do
 208
 209   do j = 2, ipph
 210     do ik = 1, idl1
 211       ch2(1,ik,1) = ch2(1,ik,1)+c2(1,ik,j)
 212     end do
 213   end do
 214
 215   if ( ido < l1 ) then
 216
 217     do i = 1, ido
 218       do k = 1, l1
 219         cc(1,i,1,k) = ch(1,i,k,1)
 220       end do
 221     end do
 222
 223   else
 224
 225     do k = 1, l1
 226       do i = 1, ido
 227         cc(1,i,1,k) = ch(1,i,k,1)
 228       end do
 229     end do
 230
 231   end if
 232
 233   do j = 2, ipph
 234     jc = ipp2 - j
 235     j2 = j+j
 236     do k = 1, l1
 237       cc(1,ido,j2-2,k) = ch(1,1,k,j)
 238       cc(1,1,j2-1,k) = ch(1,1,k,jc)
 239     end do
 240   end do
 241
 242   if ( ido == 1 ) then
 243     return
 244   end if
 245
 246   if ( nbd < l1 ) then
 247
 248     do j = 2, ipph
 249       jc = ipp2 - j
 250       j2 = j + j
 251       do i = 3, ido, 2
 252         ic = idp2 - i
 253         do k = 1, l1
 254           cc(1,i-1,j2-1,k) = ch(1,i-1,k,j)+ch(1,i-1,k,jc)
 255           cc(1,ic-1,j2-2,k) = ch(1,i-1,k,j)-ch(1,i-1,k,jc)
 256           cc(1,i,j2-1,k) = ch(1,i,k,j)+ch(1,i,k,jc)
 257           cc(1,ic,j2-2,k) = ch(1,i,k,jc)-ch(1,i,k,j)
 258         end do
 259       end do
 260    end do
 261
 262   else
 263
 264     do j = 2, ipph
 265       jc = ipp2 - j
 266       j2 = j + j
 267       do k = 1, l1
 268         do i = 3, ido, 2
 269           ic = idp2 - i
 270           cc(1,i-1,j2-1,k)  = ch(1,i-1,k,j) + ch(1,i-1,k,jc)
 271           cc(1,ic-1,j2-2,k) = ch(1,i-1,k,j) - ch(1,i-1,k,jc)
 272           cc(1,i,j2-1,k)    = ch(1,i,k,j)   + ch(1,i,k,jc)
 273           cc(1,ic,j2-2,k)   = ch(1,i,k,jc)  - ch(1,i,k,j)
 274         end do
 275       end do
 276     end do
 277
 278   end if
 279
 280   return
 281 end