fftsg3d.f

2维fft程序

              end do
          else if (n1 .eq. 4) then
              do j = 0, n2 - 1
                  t(2 * j) = a(0, j, k)
                  t(2 * j + 1) = a(1, j, k)
                  t(2 * n2 + 2 * j) = a(2, j, k)
                  t(2 * n2 + 2 * j + 1) = a(3, j, k)
              end do
              call cdft(2 * n2, isgn, t, ip, w)
              call cdft(2 * n2, isgn, t(2 * n2), ip, w)
              do j = 0, n2 - 1
                  a(0, j, k) = t(2 * j)
                  a(1, j, k) = t(2 * j + 1)
                  a(2, j, k) = t(2 * n2 + 2 * j)
                  a(3, j, k) = t(2 * n2 + 2 * j + 1)
              end do
          else if (n1 .eq. 2) then
              do j = 0, n2 - 1
                  t(2 * j) = a(0, j, k)
                  t(2 * j + 1) = a(1, j, k)
              end do
              call cdft(2 * n2, isgn, t, ip, w)
              do j = 0, n2 - 1
                  a(0, j, k) = t(2 * j)
                  a(1, j, k) = t(2 * j + 1)
              end do
          end if
          if (icr .ne. 0 .and. isgn .lt. 0) then
              do j = 0, n2 - 1
                  call rdft(n1, isgn, a(0, j, k), ip, w)
              end do
          end if
      end do
      end
!
      subroutine cdft3db_sub(n1max, n2max, n1, n2, n3, 
     &    isgn, a, t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 
     &    i, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      if (n1 .gt. 4) then
          do j = 0, n2 - 1
              do i = 0, n1 - 8, 8
                  do k = 0, n3 - 1
                      t(2 * k) = a(i, j, k)
                      t(2 * k + 1) = a(i + 1, j, k)
                      t(2 * n3 + 2 * k) = a(i + 2, j, k)
                      t(2 * n3 + 2 * k + 1) = a(i + 3, j, k)
                      t(4 * n3 + 2 * k) = a(i + 4, j, k)
                      t(4 * n3 + 2 * k + 1) = a(i + 5, j, k)
                      t(6 * n3 + 2 * k) = a(i + 6, j, k)
                      t(6 * n3 + 2 * k + 1) = a(i + 7, j, k)
                  end do
                  call cdft(2 * n3, isgn, t, ip, w)
                  call cdft(2 * n3, isgn, t(2 * n3), ip, w)
                  call cdft(2 * n3, isgn, t(4 * n3), ip, w)
                  call cdft(2 * n3, isgn, t(6 * n3), ip, w)
                  do k = 0, n3 - 1
                      a(i, j, k) = t(2 * k)
                      a(i + 1, j, k) = t(2 * k + 1)
                      a(i + 2, j, k) = t(2 * n3 + 2 * k)
                      a(i + 3, j, k) = t(2 * n3 + 2 * k + 1)
                      a(i + 4, j, k) = t(4 * n3 + 2 * k)
                      a(i + 5, j, k) = t(4 * n3 + 2 * k + 1)
                      a(i + 6, j, k) = t(6 * n3 + 2 * k)
                      a(i + 7, j, k) = t(6 * n3 + 2 * k + 1)
                  end do
              end do
          end do
      else if (n1 .eq. 4) then
          do j = 0, n2 - 1
              do k = 0, n3 - 1
                  t(2 * k) = a(0, j, k)
                  t(2 * k + 1) = a(1, j, k)
                  t(2 * n3 + 2 * k) = a(2, j, k)
                  t(2 * n3 + 2 * k + 1) = a(3, j, k)
              end do
              call cdft(2 * n3, isgn, t, ip, w)
              call cdft(2 * n3, isgn, t(2 * n3), ip, w)
              do k = 0, n3 - 1
                  a(0, j, k) = t(2 * k)
                  a(1, j, k) = t(2 * k + 1)
                  a(2, j, k) = t(2 * n3 + 2 * k)
                  a(3, j, k) = t(2 * n3 + 2 * k + 1)
              end do
          end do
      else if (n1 .eq. 2) then
          do j = 0, n2 - 1
              do k = 0, n3 - 1
                  t(2 * k) = a(0, j, k)
                  t(2 * k + 1) = a(1, j, k)
              end do
              call cdft(2 * n3, isgn, t, ip, w)
              do k = 0, n3 - 1
                  a(0, j, k) = t(2 * k)
                  a(1, j, k) = t(2 * k + 1)
              end do
          end do
      end if
      end
!
      subroutine rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a)
      integer n1max, n2max, n1, n2, n3, isgn, 
     &    n2h, n3h, i, j, k, l
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), xi
      n2h = n2 / 2
      n3h = n3 / 2
      if (isgn .lt. 0) then
          do k = 1, n3h - 1
              l = n3 - k
              xi = a(0, 0, k) - a(0, 0, l)
              a(0, 0, k) = a(0, 0, k) + a(0, 0, l)
              a(0, 0, l) = xi
              xi = a(1, 0, l) - a(1, 0, k)
              a(1, 0, k) = a(1, 0, k) + a(1, 0, l)
              a(1, 0, l) = xi
              xi = a(0, n2h, k) - a(0, n2h, l)
              a(0, n2h, k) = a(0, n2h, k) + a(0, n2h, l)
              a(0, n2h, l) = xi
              xi = a(1, n2h, l) - a(1, n2h, k)
              a(1, n2h, k) = a(1, n2h, k) + a(1, n2h, l)
              a(1, n2h, l) = xi
              do i = 1, n2h - 1
                  j = n2 - i
                  xi = a(0, i, k) - a(0, j, l)
                  a(0, i, k) = a(0, i, k) + a(0, j, l)
                  a(0, j, l) = xi
                  xi = a(1, j, l) - a(1, i, k)
                  a(1, i, k) = a(1, i, k) + a(1, j, l)
                  a(1, j, l) = xi
                  xi = a(0, i, l) - a(0, j, k)
                  a(0, i, l) = a(0, i, l) + a(0, j, k)
                  a(0, j, k) = xi
                  xi = a(1, j, k) - a(1, i, l)
                  a(1, i, l) = a(1, i, l) + a(1, j, k)
                  a(1, j, k) = xi
              end do
          end do
          do i = 1, n2h - 1
              j = n2 - i
              xi = a(0, i, 0) - a(0, j, 0)
              a(0, i, 0) = a(0, i, 0) + a(0, j, 0)
              a(0, j, 0) = xi
              xi = a(1, j, 0) - a(1, i, 0)
              a(1, i, 0) = a(1, i, 0) + a(1, j, 0)
              a(1, j, 0) = xi
              xi = a(0, i, n3h) - a(0, j, n3h)
              a(0, i, n3h) = a(0, i, n3h) + a(0, j, n3h)
              a(0, j, n3h) = xi
              xi = a(1, j, n3h) - a(1, i, n3h)
              a(1, i, n3h) = a(1, i, n3h) + a(1, j, n3h)
              a(1, j, n3h) = xi
          end do
      else
          do k = 1, n3h - 1
              l = n3 - k
              a(0, 0, l) = 0.5d0 * (a(0, 0, k) - a(0, 0, l))
              a(0, 0, k) = a(0, 0, k) - a(0, 0, l)
              a(1, 0, l) = 0.5d0 * (a(1, 0, k) + a(1, 0, l))
              a(1, 0, k) = a(1, 0, k) - a(1, 0, l)
              a(0, n2h, l) = 0.5d0 * (a(0, n2h, k) - a(0, n2h, l))
              a(0, n2h, k) = a(0, n2h, k) - a(0, n2h, l)
              a(1, n2h, l) = 0.5d0 * (a(1, n2h, k) + a(1, n2h, l))
              a(1, n2h, k) = a(1, n2h, k) - a(1, n2h, l)
              do i = 1, n2h - 1
                  j = n2 - i
                  a(0, j, l) = 0.5d0 * (a(0, i, k) - a(0, j, l))
                  a(0, i, k) = a(0, i, k) - a(0, j, l)
                  a(1, j, l) = 0.5d0 * (a(1, i, k) + a(1, j, l))
                  a(1, i, k) = a(1, i, k) - a(1, j, l)
                  a(0, j, k) = 0.5d0 * (a(0, i, l) - a(0, j, k))
                  a(0, i, l) = a(0, i, l) - a(0, j, k)
                  a(1, j, k) = 0.5d0 * (a(1, i, l) + a(1, j, k))
                  a(1, i, l) = a(1, i, l) - a(1, j, k)
              end do
          end do
          do i = 1, n2h - 1
              j = n2 - i
              a(0, j, 0) = 0.5d0 * (a(0, i, 0) - a(0, j, 0))
              a(0, i, 0) = a(0, i, 0) - a(0, j, 0)
              a(1, j, 0) = 0.5d0 * (a(1, i, 0) + a(1, j, 0))
              a(1, i, 0) = a(1, i, 0) - a(1, j, 0)
              a(0, j, n3h) = 0.5d0 * (a(0, i, n3h) - a(0, j, n3h))
              a(0, i, n3h) = a(0, i, n3h) - a(0, j, n3h)
              a(1, j, n3h) = 0.5d0 * (a(1, i, n3h) + a(1, j, n3h))
              a(1, i, n3h) = a(1, i, n3h) - a(1, j, n3h)
          end do
      end if
      end
!
      subroutine ddxt3da_sub(n1max, n2max, n1, n2, n3, ics, 
     &    isgn, a, t, ip, w)
      integer n1max, n2max, n1, n2, n3, ics, 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 (ics .eq. 0) then
              do j = 0, n2 - 1
                  call ddct(n1, isgn, a(0, j, k), ip, w)
              end do
          else
              do j = 0, n2 - 1
                  call ddst(n1, isgn, a(0, j, k), ip, w)
              end do
          end if
          if (n1 .gt. 2) then
              do i = 0, n1 - 4, 4
                  do j = 0, n2 - 1
                      t(j) = a(i, j, k)
                      t(n2 + j) = a(i + 1, j, k)
                      t(2 * n2 + j) = a(i + 2, j, k)
                      t(3 * n2 + j) = a(i + 3, j, k)
                  end do
                  if (ics .eq. 0) then
                      call ddct(n2, isgn, t, ip, w)
                      call ddct(n2, isgn, t(n2), ip, w)
                      call ddct(n2, isgn, t(2 * n2), ip, w)
                      call ddct(n2, isgn, t(3 * n2), ip, w)
                  else
                      call ddst(n2, isgn, t, ip, w)
                      call ddst(n2, isgn, t(n2), ip, w)
                      call ddst(n2, isgn, t(2 * n2), ip, w)
                      call ddst(n2, isgn, t(3 * n2), ip, w)
                  end if
                  do j = 0, n2 - 1
                      a(i, j, k) = t(j)
                      a(i + 1, j, k) = t(n2 + j)
                      a(i + 2, j, k) = t(2 * n2 + j)
                      a(i + 3, j, k) = t(3 * n2 + j)
                  end do
              end do
          else if (n1 .eq. 2) then
              do j = 0, n2 - 1
                  t(j) = a(0, j, k)
                  t(n2 + j) = a(1, j, k)
              end do
              if (ics .eq. 0) then
                  call ddct(n2, isgn, t, ip, w)
                  call ddct(n2, isgn, t(n2), ip, w)
              else
                  call ddst(n2, isgn, t, ip, w)
                  call ddst(n2, isgn, t(n2), ip, w)
              end if
              do j = 0, n2 - 1
                  a(0, j, k) = t(j)
                  a(1, j, k) = t(n2 + j)
              end do
          end if
      end do
      end
!
      subroutine ddxt3db_sub(n1max, n2max, n1, n2, n3, ics, 
     &    isgn, a, t, ip, w)
      integer n1max, n2max, n1, n2, n3, ics, isgn, 
     &    ip(0 : *), i, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      if (n1 .gt. 2) then
          do j = 0, n2 - 1
              do i = 0, n1 - 4, 4
                  do k = 0, n3 - 1
                      t(k) = a(i, j, k)
                      t(n3 + k) = a(i + 1, j, k)
                      t(2 * n3 + k) = a(i + 2, j, k)
                      t(3 * n3 + k) = a(i + 3, j, k)
                  end do
                  if (ics .eq. 0) then
                      call ddct(n3, isgn, t, ip, w)
                      call ddct(n3, isgn, t(n3), ip, w)
                      call ddct(n3, isgn, t(2 * n3), ip, w)
                      call ddct(n3, isgn, t(3 * n3), ip, w)
                  else
                      call ddst(n3, isgn, t, ip, w)
                      call ddst(n3, isgn, t(n3), ip, w)
                      call ddst(n3, isgn, t(2 * n3), ip, w)
                      call ddst(n3, isgn, t(3 * n3), ip, w)
                  end if
                  do k = 0, n3 - 1
                      a(i, j, k) = t(k)
                      a(i + 1, j, k) = t(n3 + k)
                      a(i + 2, j, k) = t(2 * n3 + k)
                      a(i + 3, j, k) = t(3 * n3 + k)
                  end do
              end do
          end do
      else if (n1 .eq. 2) then
          do j = 0, n2 - 1
              do k = 0, n3 - 1
                  t(k) = a(0, j, k)
                  t(n3 + k) = a(1, j, k)
              end do
              if (ics .eq. 0) then
                  call ddct(n3, isgn, t, ip, w)
                  call ddct(n3, isgn, t(n3), ip, w)
              else
                  call ddst(n3, isgn, t, ip, w)
                  call ddst(n3, isgn, t(n3), ip, w)
              end if
              do k = 0, n3 - 1
                  a(0, j, k) = t(k)
                  a(1, j, k) = t(n3 + k)
              end do
          end do
      end if
      end
!