external/fftpack/fftpack5/c1fgkb.F

   1 subroutine c1fgkb ( ido, ip, l1, lid, na, cc, cc1, in1, ch, ch1, in2, wa )
   2
   3 !*****************************************************************************80
   4 !
   5 !! C1FGKB 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 ) ido
  37   integer ( kind = 4 ) in1
  38   integer ( kind = 4 ) in2
  39   integer ( kind = 4 ) ip
  40   integer ( kind = 4 ) l1
  41   integer ( kind = 4 ) lid
  42
  43   real ( kind = 4 ) cc(in1,l1,ip,ido)
  44   real ( kind = 4 ) cc1(in1,lid,ip)
  45   real ( kind = 4 ) ch(in2,l1,ido,ip)
  46   real ( kind = 4 ) ch1(in2,lid,ip)
  47   real ( kind = 4 ) chold1
  48   real ( kind = 4 ) chold2
  49   integer ( kind = 4 ) i
  50   integer ( kind = 4 ) idlj
  51   integer ( kind = 4 ) ipp2
  52   integer ( kind = 4 ) ipph
  53   integer ( kind = 4 ) j
  54   integer ( kind = 4 ) jc
  55   integer ( kind = 4 ) k
  56   integer ( kind = 4 ) ki
  57   integer ( kind = 4 ) l
  58   integer ( kind = 4 ) lc
  59   integer ( kind = 4 ) na
  60   real ( kind = 4 ) wa(ido,ip-1,2)
  61   real ( kind = 4 ) wai
  62   real ( kind = 4 ) war
  63
  64   ipp2 = ip + 2
  65   ipph = ( ip + 1 ) / 2
  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      do ki = 1, lid
  92        cc1(1,ki,l) = ch1(1,ki,1)+wa(1,l-1,1)*ch1(1,ki,2)
  93        cc1(1,ki,lc) = wa(1,l-1,2)*ch1(1,ki,ip)
  94        cc1(2,ki,l) = ch1(2,ki,1)+wa(1,l-1,1)*ch1(2,ki,2)
  95        cc1(2,ki,lc) = wa(1,l-1,2)*ch1(2,ki,ip)
  96      end do
  97
  98      do j = 3, ipph
  99        jc = ipp2 - j
 100        idlj = mod ( ( l - 1 ) * ( j - 1 ), ip )
 101        war = wa(1,idlj,1)
 102        wai = wa(1,idlj,2)
 103        do ki = 1, lid
 104          cc1(1,ki,l) = cc1(1,ki,l)+war*ch1(1,ki,j)
 105          cc1(1,ki,lc) = cc1(1,ki,lc)+wai*ch1(1,ki,jc)
 106          cc1(2,ki,l) = cc1(2,ki,l)+war*ch1(2,ki,j)
 107          cc1(2,ki,lc) = cc1(2,ki,lc)+wai*ch1(2,ki,jc)
 108        end do
 109      end do
 110
 111   end do
 112
 113   if ( 1 < ido .or. na == 1 ) then
 114
 115     do ki = 1, lid
 116       ch1(1,ki,1) = cc1(1,ki,1)
 117       ch1(2,ki,1) = cc1(2,ki,1)
 118     end do
 119
 120     do j = 2, ipph
 121       jc = ipp2 - j
 122       do ki = 1, lid
 123         ch1(1,ki,j) = cc1(1,ki,j)-cc1(2,ki,jc)
 124         ch1(1,ki,jc) = cc1(1,ki,j)+cc1(2,ki,jc)
 125         ch1(2,ki,jc) = cc1(2,ki,j)-cc1(1,ki,jc)
 126         ch1(2,ki,j) = cc1(2,ki,j)+cc1(1,ki,jc)
 127       end do
 128     end do
 129
 130     if ( ido == 1 ) then
 131       return
 132     end if
 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) &
 152                        -wa(i,j-1,2)*ch(2,k,i,j)
 153           cc(2,k,j,i) = wa(i,j-1,1)*ch(2,k,i,j) &
 154                        +wa(i,j-1,2)*ch(1,k,i,j)
 155         end do
 156       end do
 157     end do
 158
 159   else
 160
 161     do j = 2, ipph
 162       jc = ipp2 - j
 163       do ki = 1, lid
 164         chold1 = cc1(1,ki,j)-cc1(2,ki,jc)
 165         chold2 = cc1(1,ki,j)+cc1(2,ki,jc)
 166         cc1(1,ki,j) = chold1
 167         cc1(2,ki,jc) = cc1(2,ki,j)-cc1(1,ki,jc)
 168         cc1(2,ki,j) = cc1(2,ki,j)+cc1(1,ki,jc)
 169         cc1(1,ki,jc) = chold2
 170       end do
 171     end do
 172
 173   end if
 174
 175   return
 176 end