external/fftpack/fftpack5/d1fgkb.F

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