external/fftpack/fftpack5/z1fgkf.F

   1 subroutine z1fgkf ( ido, ip, l1, lid, na, cc, cc1, in1, ch, ch1, in2, wa )
   2
   3 !*****************************************************************************80
   4 !
   5 !! Z1FGKF is an FFTPACK5 auxiliary routine.
   6 !
   7 !
   8 !
   9 !  Modified:
  10 !
  11 !    26 Ausust 2009
  12 !
  13 !  Author:
  14 !
  15 !    Original complex single precision by Paul Swarztrauber, Richard Valent.
  16 !    Complex 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 ) ido
  35   integer ( kind = 4 ) in1
  36   integer ( kind = 4 ) in2
  37   integer ( kind = 4 ) ip
  38   integer ( kind = 4 ) l1
  39   integer ( kind = 4 ) lid
  40
  41   real ( kind = 8 ) cc(in1,l1,ip,ido)
  42   real ( kind = 8 ) cc1(in1,lid,ip)
  43   real ( kind = 8 ) ch(in2,l1,ido,ip)
  44   real ( kind = 8 ) ch1(in2,lid,ip)
  45   real ( kind = 8 ) chold1
  46   real ( kind = 8 ) chold2
  47   integer ( kind = 4 ) i
  48   integer ( kind = 4 ) idlj
  49   integer ( kind = 4 ) ipp2
  50   integer ( kind = 4 ) ipph
  51   integer ( kind = 4 ) j
  52   integer ( kind = 4 ) jc
  53   integer ( kind = 4 ) k
  54   integer ( kind = 4 ) ki
  55   integer ( kind = 4 ) l
  56   integer ( kind = 4 ) lc
  57   integer ( kind = 4 ) na
  58   real ( kind = 8 ) sn
  59   real ( kind = 8 ) wa(ido,ip-1,2)
  60   real ( kind = 8 ) wai
  61   real ( kind = 8 ) war
  62
  63   ipp2 = ip+2
  64   ipph = (ip+1)/2
  65
  66   do ki = 1, lid
  67     ch1(1,ki,1) = cc1(1,ki,1)
  68     ch1(2,ki,1) = cc1(2,ki,1)
  69   end do
  70
  71   do j = 2, ipph
  72     jc = ipp2 - j
  73     do ki = 1, lid
  74       ch1(1,ki,j) =  cc1(1,ki,j)+cc1(1,ki,jc)
  75       ch1(1,ki,jc) = cc1(1,ki,j)-cc1(1,ki,jc)
  76       ch1(2,ki,j) =  cc1(2,ki,j)+cc1(2,ki,jc)
  77       ch1(2,ki,jc) = cc1(2,ki,j)-cc1(2,ki,jc)
  78     end do
  79   end do
  80
  81   do j = 2, ipph
  82     do ki = 1, lid
  83       cc1(1,ki,1) = cc1(1,ki,1) + ch1(1,ki,j)
  84       cc1(2,ki,1) = cc1(2,ki,1) + ch1(2,ki,j)
  85     end do
  86   end do
  87
  88   do l = 2, ipph
  89
  90     lc = ipp2 - l
  91
  92     do ki = 1, lid
  93       cc1(1,ki,l)  = ch1(1,ki,1) + wa(1,l-1,1) * ch1(1,ki,2)
  94       cc1(1,ki,lc) =             - wa(1,l-1,2) * ch1(1,ki,ip)
  95       cc1(2,ki,l)  = ch1(2,ki,1) + wa(1,l-1,1) * ch1(2,ki,2)
  96       cc1(2,ki,lc) =             - wa(1,l-1,2) * ch1(2,ki,ip)
  97     end do
  98
  99     do j = 3, ipph
 100
 101       jc = ipp2 - j
 102       idlj = mod ( ( l - 1 ) * ( j - 1 ), ip )
 103       war = wa(1,idlj,1)
 104       wai = -wa(1,idlj,2)
 105
 106       do ki = 1, lid
 107         cc1(1,ki,l) = cc1(1,ki,l)+war*ch1(1,ki,j)
 108         cc1(1,ki,lc) = cc1(1,ki,lc)+wai*ch1(1,ki,jc)
 109         cc1(2,ki,l) = cc1(2,ki,l)+war*ch1(2,ki,j)
 110         cc1(2,ki,lc) = cc1(2,ki,lc)+wai*ch1(2,ki,jc)
 111       end do
 112
 113     end do
 114
 115   end do
 116
 117   if ( 1 < ido ) then
 118
 119     do ki = 1, lid
 120       ch1(1,ki,1) = cc1(1,ki,1)
 121       ch1(2,ki,1) = cc1(2,ki,1)
 122     end do
 123
 124     do j = 2, ipph
 125       jc = ipp2 - j
 126       do ki = 1, lid
 127         ch1(1,ki,j) = cc1(1,ki,j)-cc1(2,ki,jc)
 128         ch1(2,ki,j) = cc1(2,ki,j)+cc1(1,ki,jc)
 129         ch1(1,ki,jc) = cc1(1,ki,j)+cc1(2,ki,jc)
 130         ch1(2,ki,jc) = cc1(2,ki,j)-cc1(1,ki,jc)
 131       end do
 132     end do
 133
 134     do i = 1, ido
 135       do k = 1, l1
 136         cc(1,k,1,i) = ch(1,k,i,1)
 137         cc(2,k,1,i) = ch(2,k,i,1)
 138       end do
 139     end do
 140
 141     do j = 2, ip
 142       do k = 1, l1
 143         cc(1,k,j,1) = ch(1,k,1,j)
 144         cc(2,k,j,1) = ch(2,k,1,j)
 145       end do
 146     end do
 147
 148     do j = 2, ip
 149       do i = 2, ido
 150         do k = 1, l1
 151           cc(1,k,j,i) = wa(i,j-1,1)*ch(1,k,i,j) + wa(i,j-1,2)*ch(2,k,i,j)
 152           cc(2,k,j,i) = wa(i,j-1,1)*ch(2,k,i,j) - wa(i,j-1,2)*ch(1,k,i,j)
 153         end do
 154       end do
 155     end do
 156
 157   else if ( na == 1 ) then
 158
 159     sn = 1.0D+00 / real ( ip * l1, kind = 8 )
 160
 161     do ki = 1, lid
 162       ch1(1,ki,1) = sn * cc1(1,ki,1)
 163       ch1(2,ki,1) = sn * cc1(2,ki,1)
 164     end do
 165
 166     do j = 2, ipph
 167       jc = ipp2 - j
 168       do ki = 1, lid
 169         ch1(1,ki,j) = sn*(cc1(1,ki,j)-cc1(2,ki,jc))
 170         ch1(2,ki,j) = sn*(cc1(2,ki,j)+cc1(1,ki,jc))
 171         ch1(1,ki,jc) = sn*(cc1(1,ki,j)+cc1(2,ki,jc))
 172         ch1(2,ki,jc) = sn*(cc1(2,ki,j)-cc1(1,ki,jc))
 173       end do
 174     end do
 175
 176   else
 177
 178     sn = 1.0D+00 / real ( ip * l1, kind = 8 )
 179
 180     do ki = 1, lid
 181       cc1(1,ki,1) = sn * cc1(1,ki,1)
 182       cc1(2,ki,1) = sn * cc1(2,ki,1)
 183     end do
 184
 185     do j = 2, ipph
 186       jc = ipp2 - j
 187       do ki = 1, lid
 188         chold1 = sn*(cc1(1,ki,j)-cc1(2,ki,jc))
 189         chold2 = sn*(cc1(1,ki,j)+cc1(2,ki,jc))
 190         cc1(1,ki,j) = chold1
 191         cc1(2,ki,jc) = sn*(cc1(2,ki,j)-cc1(1,ki,jc))
 192         cc1(2,ki,j) = sn*(cc1(2,ki,j)+cc1(1,ki,jc))
 193         cc1(1,ki,jc) = chold2
 194       end do
 195     end do
 196
 197   end if
 198
 199   return
 200 end