external/fftpack/fftpack5/c1fgkf.F

   1 subroutine c1fgkf ( ido, ip, l1, lid, na, cc, cc1, in1, ch, ch1, in2, wa )
   2
   3 !*****************************************************************************80
   4 !
   5 !! C1FGKF 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 ) sn
  61   real ( kind = 4 ) wa(ido,ip-1,2)
  62   real ( kind = 4 ) wai
  63   real ( kind = 4 ) war
  64
  65   ipp2 = ip+2
  66   ipph = (ip+1)/2
  67
  68   do ki = 1, lid
  69     ch1(1,ki,1) = cc1(1,ki,1)
  70     ch1(2,ki,1) = cc1(2,ki,1)
  71   end do
  72
  73   do j = 2, ipph
  74     jc = ipp2 - j
  75     do ki = 1, lid
  76       ch1(1,ki,j) =  cc1(1,ki,j)+cc1(1,ki,jc)
  77       ch1(1,ki,jc) = cc1(1,ki,j)-cc1(1,ki,jc)
  78       ch1(2,ki,j) =  cc1(2,ki,j)+cc1(2,ki,jc)
  79       ch1(2,ki,jc) = cc1(2,ki,j)-cc1(2,ki,jc)
  80     end do
  81   end do
  82
  83   do j = 2, ipph
  84     do ki = 1, lid
  85       cc1(1,ki,1) = cc1(1,ki,1) + ch1(1,ki,j)
  86       cc1(2,ki,1) = cc1(2,ki,1) + ch1(2,ki,j)
  87     end do
  88   end do
  89
  90   do l = 2, ipph
  91
  92     lc = ipp2 - l
  93
  94     do ki = 1, lid
  95       cc1(1,ki,l)  = ch1(1,ki,1) + wa(1,l-1,1) * ch1(1,ki,2)
  96       cc1(1,ki,lc) =             - wa(1,l-1,2) * ch1(1,ki,ip)
  97       cc1(2,ki,l)  = ch1(2,ki,1) + wa(1,l-1,1) * ch1(2,ki,2)
  98       cc1(2,ki,lc) =             - wa(1,l-1,2) * ch1(2,ki,ip)
  99     end do
 100
 101     do j = 3, ipph
 102
 103       jc = ipp2 - j
 104       idlj = mod ( ( l - 1 ) * ( j - 1 ), ip )
 105       war = wa(1,idlj,1)
 106       wai = -wa(1,idlj,2)
 107
 108       do ki = 1, lid
 109         cc1(1,ki,l) = cc1(1,ki,l)+war*ch1(1,ki,j)
 110         cc1(1,ki,lc) = cc1(1,ki,lc)+wai*ch1(1,ki,jc)
 111         cc1(2,ki,l) = cc1(2,ki,l)+war*ch1(2,ki,j)
 112         cc1(2,ki,lc) = cc1(2,ki,lc)+wai*ch1(2,ki,jc)
 113       end do
 114
 115     end do
 116
 117   end do
 118
 119   if ( 1 < ido ) then
 120
 121     do ki = 1, lid
 122       ch1(1,ki,1) = cc1(1,ki,1)
 123       ch1(2,ki,1) = cc1(2,ki,1)
 124     end do
 125
 126     do j = 2, ipph
 127       jc = ipp2 - j
 128       do ki = 1, lid
 129         ch1(1,ki,j) = cc1(1,ki,j)-cc1(2,ki,jc)
 130         ch1(2,ki,j) = cc1(2,ki,j)+cc1(1,ki,jc)
 131         ch1(1,ki,jc) = cc1(1,ki,j)+cc1(2,ki,jc)
 132         ch1(2,ki,jc) = cc1(2,ki,j)-cc1(1,ki,jc)
 133       end do
 134     end do
 135
 136     do i = 1, ido
 137       do k = 1, l1
 138         cc(1,k,1,i) = ch(1,k,i,1)
 139         cc(2,k,1,i) = ch(2,k,i,1)
 140       end do
 141     end do
 142
 143     do j = 2, ip
 144       do k = 1, l1
 145         cc(1,k,j,1) = ch(1,k,1,j)
 146         cc(2,k,j,1) = ch(2,k,1,j)
 147       end do
 148     end do
 149
 150     do j = 2, ip
 151       do i = 2, ido
 152         do k = 1, l1
 153           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)
 154           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)
 155         end do
 156       end do
 157     end do
 158
 159   else if ( na == 1 ) then
 160
 161     sn = 1.0E+00 / real ( ip * l1, kind = 4 )
 162
 163     do ki = 1, lid
 164       ch1(1,ki,1) = sn * cc1(1,ki,1)
 165       ch1(2,ki,1) = sn * cc1(2,ki,1)
 166     end do
 167
 168     do j = 2, ipph
 169       jc = ipp2 - j
 170       do ki = 1, lid
 171         ch1(1,ki,j) =  sn * ( cc1(1,ki,j) - cc1(2,ki,jc) )
 172         ch1(2,ki,j) =  sn * ( cc1(2,ki,j) + cc1(1,ki,jc) )
 173         ch1(1,ki,jc) = sn * ( cc1(1,ki,j) + cc1(2,ki,jc) )
 174         ch1(2,ki,jc) = sn * ( cc1(2,ki,j) - cc1(1,ki,jc) )
 175       end do
 176     end do
 177
 178   else
 179
 180     sn = 1.0E+00 / real ( ip * l1, kind = 4 )
 181
 182     do ki = 1, lid
 183       cc1(1,ki,1) = sn * cc1(1,ki,1)
 184       cc1(2,ki,1) = sn * cc1(2,ki,1)
 185     end do
 186
 187     do j = 2, ipph
 188       jc = ipp2 - j
 189       do ki = 1, lid
 190         chold1 = sn*(cc1(1,ki,j)-cc1(2,ki,jc))
 191         chold2 = sn*(cc1(1,ki,j)+cc1(2,ki,jc))
 192         cc1(1,ki,j) = chold1
 193         cc1(2,ki,jc) = sn*(cc1(2,ki,j)-cc1(1,ki,jc))
 194         cc1(2,ki,j) = sn*(cc1(2,ki,j)+cc1(1,ki,jc))
 195         cc1(1,ki,jc) = chold2
 196       end do
 197     end do
 198
 199   end if
 200
 201   return
 202 end