fftsg3d.f

2维fft程序

!                     do j1 = 0, n1 - 1
!                         a(j1,j2,j3) = a(j1,j2,j3) * (8.0d0/n1/n2/n3)
!                     end do
!                 end do
!             end do
!         .
!
!
! -------- DST (Discrete Sine Transform) / Inverse of DST --------
!     [definition]
!         <case1> IDST (excluding scale)
!             S(k1,k2,k3) = sum_j1=1^n1 sum_j2=1^n2 sum_j3=1^n3
!                               A(j1,j2,j3) * 
!                               sin(pi*j1*(k1+1/2)/n1) * 
!                               sin(pi*j2*(k2+1/2)/n2) * 
!                               sin(pi*j3*(k3+1/2)/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         <case2> DST
!             S(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               sin(pi*(j1+1/2)*k1/n1) * 
!                               sin(pi*(j2+1/2)*k2/n2) * 
!                               sin(pi*(j3+1/2)*k3/n3), 
!                               0<k1<=n1, 0<k2<=n2, 0<k3<=n3
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call ddst3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call ddst3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!     [parameters]
!         n1max  :row1 size of the 3D array (integer)
!         n2max  :row2 size of the 3D array (integer)
!         n1     :data length (integer)
!                 n1 >= 2, n1 = power of 2
!         n2     :data length (integer)
!                 n2 >= 2, n2 = power of 2
!         n3     :data length (integer)
!                 n3 >= 2, n3 = power of 2
!         a(0:n1-1,0:n2-1,0:n3-1)
!                :input/output data (real*8)
!                 <case1>
!                     input data
!                         a(mod(j1,n1),mod(j2,n2),mod(j3,n3)) = A(j1,j2,j3), 
!                             0<j1<=n1, 0<j2<=n2, 0<j3<=n3
!                     output data
!                         a(k1,k2,k3) = S(k1,k2,k3), 
!                             0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!                 <case2>
!                     output data
!                         a(mod(k1,n1),mod(k2,n2),mod(k3,n3)) = S(k1,k2,k3), 
!                             0<k1<=n1, 0<k2<=n2, 0<k3<=n3
!         t(0:*) :work area (real*8)
!                 length of t >= max(4*n2, 4*n3)
!         ip(0:*):work area for bit reversal (integer)
!                 length of ip >= 2+sqrt(n)
!                 (n = max(n1/2, n2/2, n3/2))
!                 ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:*) :cos/sin table (real*8)
!                 length of w >= max(n1*3/2, n2*3/2, n3*3/2)
!                 w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of 
!             call ddst3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!         is 
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     a(0, j2, j3) = a(0, j2, j3) * 0.5d0
!                 end do
!                 do j1 = 0, n1 - 1
!                     a(j1, 0, j3) = a(j1, 0, j3) * 0.5d0
!                 end do
!             end do
!             do j2 = 0, n2 - 1
!                 do j1 = 0, n1 - 1
!                     a(j1, j2, 0) = a(j1, j2, 0) * 0.5d0
!                 end do
!             end do
!             call ddst3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     do j1 = 0, n1 - 1
!                         a(j1,j2,j3) = a(j1,j2,j3) * (8.0d0/n1/n2/n3)
!                     end do
!                 end do
!             end do
!         .
!
!
      subroutine cdft3d(n1max, n2max, n1, n2, n3, isgn, a, 
     &    t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), n
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      n = 2 * max(n2, n3)
      n = max(n, n1)
      if (n .gt. 4 * ip(0)) then
          call makewt(n / 4, ip, w)
      end if
      call xdft3da_sub(n1max, n2max, n1, n2, n3, 0, 
     &    isgn, a, t, ip, w)
      call cdft3db_sub(n1max, n2max, n1, n2, n3, 
     &    isgn, a, t, ip, w)
      end
!
      subroutine rdft3d(n1max, n2max, n1, n2, n3, isgn, a, 
     &    t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 
     &    n, nw, nc
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      n = 2 * max(n2, n3)
      n = max(n, n1)
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n1 .gt. 4 * nc) then
          nc = n1 / 4
          call makect(nc, ip, w(nw))
      end if
      if (isgn .lt. 0) then
          call rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a)
          call cdft3db_sub(n1max, n2max, n1, n2, n3, 
     &        isgn, a, t, ip, w)
          call xdft3da_sub(n1max, n2max, n1, n2, n3, 1, 
     &        isgn, a, t, ip, w)
      else
          call xdft3da_sub(n1max, n2max, n1, n2, n3, 1, 
     &        isgn, a, t, ip, w)
          call cdft3db_sub(n1max, n2max, n1, n2, n3, 
     &        isgn, a, t, ip, w)
          call rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a)
      end if
      end
!
      subroutine rdft3dsort(n1max, n2max, n1, n2, n3, isgn, a)
      integer n1max, n2max, n1, n2, n3, isgn, n2h, n3h, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), x, y
      n2h = n2 / 2
      n3h = n3 / 2
      if (isgn .lt. 0) then
          do k = 0, n3 - 1
              do j = n2h + 1, n2 - 1
                  a(0, j, k) = a(n1 + 1, j, k)
                  a(1, j, k) = a(n1, j, k)
              end do
          end do
          do k = n3h + 1, n3 - 1
              a(0, 0, k) = a(n1 + 1, 0, k)
              a(1, 0, k) = a(n1, 0, k)
              a(0, n2h, k) = a(n1 + 1, n2h, k)
              a(1, n2h, k) = a(n1, n2h, k)
          end do
          a(1, 0, 0) = a(n1, 0, 0)
          a(1, n2h, 0) = a(n1, n2h, 0)
          a(1, 0, n3h) = a(n1, 0, n3h)
          a(1, n2h, n3h) = a(n1, n2h, n3h)
      else
          do j = n2h + 1, n2 - 1
              y = a(0, j, 0)
              x = a(1, j, 0)
              a(n1, j, 0) = x
              a(n1 + 1, j, 0) = y
              a(n1, n2 - j, 0) = x
              a(n1 + 1, n2 - j, 0) = -y
              a(0, j, 0) = a(0, n2 - j, 0)
              a(1, j, 0) = -a(1, n2 - j, 0)
          end do
          do k = 1, n3 - 1
              do j = n2h + 1, n2 - 1
                  y = a(0, j, k)
                  x = a(1, j, k)
                  a(n1, j, k) = x
                  a(n1 + 1, j, k) = y
                  a(n1, n2 - j, n3 - k) = x
                  a(n1 + 1, n2 - j, n3 - k) = -y
                  a(0, j, k) = a(0, n2 - j, n3 - k)
                  a(1, j, k) = -a(1, n2 - j, n3 - k)
              end do
          end do
          do k = n3h + 1, n3 - 1
              y = a(0, 0, k)
              x = a(1, 0, k)
              a(n1, 0, k) = x
              a(n1 + 1, 0, k) = y
              a(n1, 0, n3 - k) = x
              a(n1 + 1, 0, n3 - k) = -y
              a(0, 0, k) = a(0, 0, n3 - k)
              a(1, 0, k) = -a(1, 0, n3 - k)
              y = a(0, n2h, k)
              x = a(1, n2h, k)
              a(n1, n2h, k) = x
              a(n1 + 1, n2h, k) = y
              a(n1, n2h, n3 - k) = x
              a(n1 + 1, n2h, n3 - k) = -y
              a(0, n2h, k) = a(0, n2h, n3 - k)
              a(1, n2h, k) = -a(1, n2h, n3 - k)
          end do
          a(n1, 0, 0) = a(1, 0, 0)
          a(n1 + 1, 0, 0) = 0
          a(1, 0, 0) = 0
          a(n1, n2h, 0) = a(1, n2h, 0)
          a(n1 + 1, n2h, 0) = 0
          a(1, n2h, 0) = 0
          a(n1, 0, n3h) = a(1, 0, n3h)
          a(n1 + 1, 0, n3h) = 0
          a(1, 0, n3h) = 0
          a(n1, n2h, n3h) = a(1, n2h, n3h)
          a(n1 + 1, n2h, n3h) = 0
          a(1, n2h, n3h) = 0
      end if
      end
!
      subroutine ddct3d(n1max, n2max, n1, n2, n3, isgn, a, 
     &    t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 
     &    n, nw, nc
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      n = max(n2, n3)
      n = max(n, n1)
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. nc) then
          nc = n
          call makect(nc, ip, w(nw))
      end if
      call ddxt3da_sub(n1max, n2max, n1, n2, n3, 0, 
     &    isgn, a, t, ip, w)
      call ddxt3db_sub(n1max, n2max, n1, n2, n3, 0, 
     &    isgn, a, t, ip, w)
      end
!
      subroutine ddst3d(n1max, n2max, n1, n2, n3, isgn, a, 
     &    t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 
     &    n, nw, nc
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      n = max(n2, n3)
      n = max(n, n1)
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. nc) then
          nc = n
          call makect(nc, ip, w(nw))
      end if
      call ddxt3da_sub(n1max, n2max, n1, n2, n3, 1, 
     &    isgn, a, t, ip, w)
      call ddxt3db_sub(n1max, n2max, n1, n2, n3, 1, 
     &    isgn, a, t, ip, w)
      end
!
! -------- child routines --------
!
      subroutine xdft3da_sub(n1max, n2max, n1, n2, n3, icr, 
     &    isgn, a, t, ip, w)
      integer n1max, n2max, n1, n2, n3, icr, isgn, 
     &    ip(0 : *), i, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      do k = 0, n3 - 1
          if (icr .eq. 0) then
              do j = 0, n2 - 1
                  call cdft(n1, isgn, a(0, j, k), ip, w)
              end do
          else if (isgn .ge. 0) then
              do j = 0, n2 - 1
                  call rdft(n1, isgn, a(0, j, k), ip, w)
              end do
          end if
          if (n1 .gt. 4) then
              do i = 0, n1 - 8, 8
                  do j = 0, n2 - 1
                      t(2 * j) = a(i, j, k)
                      t(2 * j + 1) = a(i + 1, j, k)
                      t(2 * n2 + 2 * j) = a(i + 2, j, k)
                      t(2 * n2 + 2 * j + 1) = a(i + 3, j, k)
                      t(4 * n2 + 2 * j) = a(i + 4, j, k)
                      t(4 * n2 + 2 * j + 1) = a(i + 5, j, k)
                      t(6 * n2 + 2 * j) = a(i + 6, j, k)
                      t(6 * n2 + 2 * j + 1) = a(i + 7, j, k)
                  end do
                  call cdft(2 * n2, isgn, t, ip, w)
                  call cdft(2 * n2, isgn, t(2 * n2), ip, w)
                  call cdft(2 * n2, isgn, t(4 * n2), ip, w)
                  call cdft(2 * n2, isgn, t(6 * n2), ip, w)
                  do j = 0, n2 - 1
                      a(i, j, k) = t(2 * j)
                      a(i + 1, j, k) = t(2 * j + 1)
                      a(i + 2, j, k) = t(2 * n2 + 2 * j)
                      a(i + 3, j, k) = t(2 * n2 + 2 * j + 1)
                      a(i + 4, j, k) = t(4 * n2 + 2 * j)
                      a(i + 5, j, k) = t(4 * n2 + 2 * j + 1)
                      a(i + 6, j, k) = t(6 * n2 + 2 * j)
                      a(i + 7, j, k) = t(6 * n2 + 2 * j + 1)
                  end do