cfft4.f

网络带宽测试工具

      Subroutine cfft4( is, m, ur, ui, xr, xi, wr, wi )
! ----------------------------------------------------------------------
! --- CFFT4 performs an N-point complex-to-complex FFT, where N = 2**M.
!     On input XR holds the Real part of the complex array;
!     XI holds the Imaginary part.
!     On output XR contains the Real part of the transformed array.
!     XI then contains the Imaginary part of the result.
!     WR and WI are work arrays that store Real and Imaginary parts
!     of intermediate results.
!     UR and UI contain rotation factors for various stages of the
!     transformation.
! 
! --- Before doing the actual transformation UR and UI must be
!     initialised by calling CFFT4 with IS = 0, and M = MX, where MX
!     is the maximum value of M that will be needed for any call of
!     CFFT4. 
!
! --- A call to CFFT4 with IS = 1 (or -1) indicates a call to perform
!     a FFT with positive (or negative) exponentials.
!
!     For any call to CFFT4 M must be M >= 2.
!     Dimensions of the arrays XR, XI, UR, UI, WR, WI are:
!
!     XR(N), XI(N), WR(N), WI(N).
!     UR(2**MX+1),       UI(2**MX+1)       for MX even.
!     UR(2*2**(MX/3)+1), UI(2*2**(MX/3)+1) for MX odd.     
!
!     In this routine we use 8-byte Real elements.
! ----------------------------------------------------------------------
      Use                    numerics
      Real(l_)            :: ur(*), ui(*), xr(*), xi(*), wr(*), wi(*)
      Integer             :: is, m

      Integer             :: i,  it, j,  j0, j1, j2, j3, j4, j5, k , kn,
     &                       ku, k0, k1, k2, k3, ln, l1, l2, l3, l4, mx,
     &                       n , nu, n1, n2, n3
      Real(l_), Parameter :: pi = 3.141592653589793_l_
      Real(l_)            :: c0r, c0i, c1r, c1i, c2r, c2i, c3r, c3i,
     &                       d0r, d0i, d1r, d1i, d2r, d2i, d3r, d3i,
     &                       t  , ti , u1r, u1i, u2r, u2i, u3r, u3i,
     &                       x1 , x2 , w1 , w2
! ----------------------------------------------------------------------
!
      If ( is == 0 ) Then
! ----------------------------------------------------------------------
! --- Initialise the U array with sines and cosines in a manner that
!     permits stride one access at each FFT iteration.

         nu    = 2**((m+1)/2*2)/3
         ur(1) = 64*nu + m
         ku    = 2
         ln    = 1

         Do j = 0, (m+1)/2 - 1
            t = pi/(2*ln)
cdir$ ivdep
            Do i = 0, ln - 1
               ti       = i*t
               ur(i+ku) = Cos( ti )
               ui(i+ku) = Sin( ti )
            End Do
            ku = ku + ln
            ln = 4*ln
         End Do
         Return
      End If
! ----------------------------------------------------------------------
! --- Transformation call to CFFT4: set initial parameters.

      n  = 2**m
      k  = ur(1)
      mx = Mod( k, 64 )
! ----------------------------------------------------------------------
! --- Check for invalid input parameters.

      If ( ( is /= 1 .AND. is /= -1 ) .OR. m < 2 .OR. m > mx ) Then
         Print 1000, is, m
         Stop
      End If
! ----------------------------------------------------------------------
!
      nu = k/64
      n1 = n/4
      n2 = 2*n1
      n3 = 3*n1
      ln = n1
! ----------------------------------------------------------------------
! --- For the first two radix-4 iterations, the inner loop increments J
!     instead of K. Non-unit stride memory access is used but the longer
!     vector lengths should compensate for that.

      j0 = 1
      j1 = j0 + 1
      j2 = j0 + 2
      j3 = j0 + 3
      k0 = j0
      k1 = k0 + n1
      k2 = k0 + n2
      k3 = k0 + n3
cdir$ ivdep
      Do j = 0, ln - 1
         c0r = xr(j+k0)
         c0i = xi(j+k0)
         c1r = xr(j+k1)
         c1i = xi(j+k1)
         c2r = xr(j+k2)
         c2i = xi(j+k2)
         c3r = xr(j+k3)
         c3i = xi(j+k3)
         d0r = c0r + c2r
         d0i = c0i + c2i
         d1r = c0r - c2r
         d1i = c0i - c2i
         d2r = c1r + c3r
         d2i = c1i + c3i
         d3r = is*(-c1i + c3i)
         d3i = is*( c1r - c3r)
         j4  = 4*j
         wr(j4+j0) = d0r + d2r
         wi(j4+j0) = d0i + d2i
         wr(j4+j1) = d1r + d3r
         wi(j4+j1) = d1i + d3i
         wr(j4+j2) = d0r - d2r
         wi(j4+j2) = d0i - d2i
         wr(j4+j3) = d1r - d3r
         wi(j4+j3) = d1i - d3i
      End Do
! ----------------------------------------------------------------------
! --- If M = 2 all transformation work is done, just copy WR/WI into
!     XR/XI (label 200). 
!     If M = 3, one radix-2 iteration still has to be done (label 300).

      If ( m == 2 ) Go To 200
      If ( m == 3 ) Go To 300

! --- Else:
! ----------------------------------------------------------------------
! --- Perform the second radix-4 iteration.

      ln = ln/4
      Do k = 0, 3
         j0  = k + 1
         j1  = j0 + 4
         j2  = j0 + 8
         j3  = j0 + 12
         k0  = j0
         k1  = k0 + n1
         k2  = k0 + n2
         k3  = k0 + n3
         u1r =    ur(k+3)
         u1i = is*ui(k+3)
         u2r = u1r*u1r - u1i*u1i
         u2i = 2.0D0*u1r*u1i
         u3r = u1r*u2r - u1i*u2i
         u3i = u1r*u2i + u1i*u2r
cdir$ ivdep
         Do j = 0, ln - 1
            j4  = 4*j
            c0r = wr(j4+k0)
            c0i = wi(j4+k0)
            w1  = wr(j4+k1)
            w2  = wi(j4+k1)
            c1r = u1r*w1 - u1i*w2
            c1i = u1r*w2 + u1i*w1
            w1  = wr(j4+k2)
            w2  = wi(j4+k2)
            c2r = u2r*w1 - u2i*w2
            c2i = u2r*w2 + u2i*w1
            w1  = wr(j4+k3)
            w2  = wi(j4+k3)
            c3r = u3r*w1 - u3i*w2
            c3i = u3r*w2 + u3i*w1
            d0r = c0r + c2r
            d0i = c0i + c2i
            d1r = c0r - c2r
            d1i = c0i - c2i
            d2r = c1r + c3r
            d2i = c1i + c3i
            d3r = is*(-c1i + c3i)
            d3i = is*( c1r - c3r)
            j5  = 16*j
            xr(j5+j0) = d0r + d2r
            xi(j5+j0) = d0i + d2i
            xr(j5+j1) = d1r + d3r
            xi(j5+j1) = d1i + d3i
            xr(j5+j2) = d0r - d2r
            xi(j5+j2) = d0i - d2i
            xr(j5+j3) = d1r - d3r
            xi(j5+j3) = d1i - d3i
         End Do
      End Do
! ----------------------------------------------------------------------
! --- Set parameters for subsequent iterations.

      ku = 7
      ln = ln/4
      l1 = 16
      l2 = 32
      l3 = 48
      l4 = 64
! ----------------------------------------------------------------------
! --- For all other radix-4 iterations, the inner loops increment K. In
!     this manner, all array references have unit stride. The radix-4 
!     iterations are performed in pairs: XR/XI to WR/WI and vice 
!     versa. This eliminates the need to copy data between X and W
!     after each iteration.
! --- If IT = M - 1 then only one radix-2 iteration needs to be done
!     (label 100).

      Do it = 4, m - 1, 4

         If ( it == m - 1 ) Go To 100

         Do j = 0, ln - 1
            j0 = j*l4 + 1
            j1 = j0 + l1
            j2 = j0 + l2
            j3 = j0 + l3
            k0 = j*l1 + 1
            k1 = k0 + n1
            k2 = k0 + n2
            k3 = k0 + n3
cdir$ ivdep
            Do k = 0, l1 - 1
               u1r =    ur(k+ku)
               u1i = is*ui(k+ku)
               u2r = u1r*u1r - u1i*u1i
               u2i = 2.0D0*u1r*u1i
               u3r = u1r*u2r - u1i*u2i
               u3i = u1r*u2i + u1i*u2r
               c0r = xr(k+k0)
               c0i = xi(k+k0)
               x1  = xr(k+k1)
               x2  = xi(k+k1)
               c1r = u1r*x1 - u1i*x2
               c1i = u1r*x2 + u1i*x1
               x1  = xr(k+k2)
               x2  = xi(k+k2)
               c2r = u2r*x1 - u2i*x2
               c2i = u2r*x2 + u2i*x1
               x1  = xr(k+k3)
               x2  = xi(k+k3)
               c3r = u3r*x1 - u3i*x2
               c3i = u3r*x2 + u3i*x1
               d0r = c0r + c2r
               d0i = c0i + c2i
               d1r = c0r - c2r
               d1i = c0i - c2i
               d2r = c1r + c3r
               d2i = c1i + c3i
               d3r = is*(-c1i + c3i)
               d3i = is*( c1r - c3r)
               wr(k+j0) = d0r + d2r
               wi(k+j0) = d0i + d2i
               wr(k+j1) = d1r + d3r
               wi(k+j1) = d1i + d3i
               wr(k+j2) = d0r - d2r
               wi(k+j2) = d0i - d2i
               wr(k+j3) = d1r - d3r
               wi(k+j3) = d1i - d3i
            End Do
         End Do

         ku = ku + l1
         ln = ln/4
         l1 = l4
         l2 = 2*l1
         l3 = 3*l1
         l4 = 4*l1
! ----------------------------------------------------------------------
! --- If IT = M - 2 the computation is complete except for copying
!                   WR/WI to XR/XI (label 200).
! --- If IT = M - 3 then only one radix-2 iteration remains to be done
!                   (label 300).

         If ( it == m - 2 ) Go To 200
         If ( it == m - 3 ) Go To 300
!
         Do j = 0, ln - 1
            j0 = j*l4 + 1
            j1 = j0 + l1
            j2 = j0 + l2
            j3 = j0 + l3
            k0 = j*l1 + 1
            k1 = k0 + n1
            k2 = k0 + n2
            k3 = k0 + n3
cdir$ ivdep
            Do k = 0, l1 - 1
               u1r =    ur(k+ku)
               u1i = is*ui(k+ku)
               u2r = u1r*u1r - u1i*u1i
               u2i = 2.0D0*u1r*u1i
               u3r = u1r*u2r - u1i*u2i
               u3i = u1r*u2i + u1i*u2r
               c0r = wr(k+k0)
               c0i = wi(k+k0)
               w1  = wr(k+k1)
               w2  = wi(k+k1)
               c1r = u1r*w1 - u1i*w2
               c1i = u1r*w2 + u1i*w1
               w1  = wr(k+k2)
               w2  = wi(k+k2)
               c2r = u2r*w1 - u2i*w2
               c2i = u2r*w2 + u2i*w1
               w1  = wr(k+k3)
               w2  = wi(k+k3)
               c3r = u3r*w1 - u3i*w2
               c3i = u3r*w2 + u3i*w1
               d0r = c0r + c2r
               d0i = c0i + c2i
               d1r = c0r - c2r
               d1i = c0i - c2i
               d2r = c1r + c3r
               d2i = c1i + c3i
               d3r = is*(-c1i + c3i)
               d3i = is*( c1r - c3r)
               xr(k+j0) = d0r + d2r
               xi(k+j0) = d0i + d2i
               xr(k+j1) = d1r + d3r
               xi(k+j1) = d1i + d3i
               xr(k+j2) = d0r - d2r
               xi(k+j2) = d0i - d2i
               xr(k+j3) = d1r - d3r
               xi(k+j3) = d1i - d3i
            End Do
         End Do
         ku = ku + l1
         ln = ln/4
         l1 = l4
         l2 = 2*l1
         l3 = 3*l1
         l4 = 4*l1
      End Do
      Go To 400
! ----------------------------------------------------------------------
! --- Perform one radix-2 iteration to complete the calculation. This
!     branch is taken when the last iteration performed above is stored
!     in XR/XI.

  100 l2 = l1/2
      j0 = 1
      j1 = j0 + l1
      j2 = j0 + l2
      j3 = j1 + l2
      k0 = 1
      k1 = k0 + l1
      k2 = k0 + l2
      k3 = k1 + l2
cdir$ ivdep
      Do k = 0, l2 - 1
         u1r =    ur(2*k + ku)
         u1i = is*ui(2*k + ku)
         c0r = xr(k+k0)
         c0i = xi(k+k0)
         x1  = xr(k+k1)
         x2  = xi(k+k1)
         c1r = u1r*x1 - u1i*x2
         c1i = u1r*x2 + u1i*x1
         wr(k+j0) = c0r + c1r
         wi(k+j0) = c0i + c1i
         wr(k+j1) = c0r - c1r
         wi(k+j1) = c0i - c1i
         u2r = -is*u1i
         u2i =  is*u1r
         c0r = xr(k+k2)
         c0i = xi(k+k2)
         x1  = xr(k+k3)
         x2  = xi(k+k3)
         c1r = u2r*x1 - u2i*x2
         c1i = u2r*x2 + u2i*x1
         wr(k+j2) = c0r + c1r
         wi(k+j2) = c0i + c1i
         wr(k+j3) = c0r - c1r
         wi(k+j3) = c0i - c1i
      End Do

  200 Do k = 1, n
         xr(k) = wr(k)
         xi(k) = wi(k)
      End Do
      Go To 400
! ----------------------------------------------------------------------
! --- Perform one radix-2 iteration to complete the calculation. This
!     branch is taken if the last iteration performed above is stored
!     in W. If M = 3, the usual settings of L1 and KU do not hold, so
!     they have to be set in this case.

  300 If ( m == 3 ) Then
         l1 = 4
         ku = 3
      End If
      l2 = l1/2
      j0 = 1
      j1 = j0 + l1
      j2 = j0 + l2
      j3 = j1 + l2
      k0 = 1
      k1 = k0 + l1
      k2 = k0 + l2
      k3 = k1 + l2
cdir$ ivdep
      Do k = 0, l2 - 1
         u1r =    ur(2*k + ku)
         u1i = is*ui(2*k + ku)
         c0r = wr(k+k0)
         c0i = wi(k+k0)
         w1  = wr(k+k1)
         w2  = wi(k+k1)
         c1r = u1r*w1 - u1i*w2
         c1i = u1r*w2 + u1i*w1
         xr(k+j0) = c0r + c1r
         xi(k+j0) = c0i + c1i
         xr(k+j1) = c0r - c1r
         xi(k+j1) = c0i - c1i
         u2r = -is*u1i
         u2i =  is*u1r
         c0r = wr(k+k2)
         c0i = wi(k+k2)
         w1  = wr(k+k3)
         w2  = wi(k+k3)
         c1r = u2r*w1 - u2i*w2
         c1i = u2r*w2 + u2i*w1
         xr(k+j2) = c0r + c1r
         xi(k+j2) = c0i + c1i
         xr(k+j3) = c0r - c1r
         xi(k+j3) = c0i - c1i
      End Do
! ----------------------------------------------------------------------
! --- Finished: Return
 
  400 Return
! ----------------------------------------------------------------------
! --- Formats:

 1000 Format(' *** Error in CFFT4: Either U has not been initialised',/
     &       ,' or one of the input parameters IS or M is invalid:',/
     &       ,' IS = ', I5, '; M = ', I5)
! ----------------------------------------------------------------------
      End Subroutine cfft4