24小时热门版块排行榜

返回列表

panda0927

木虫 (初入文坛)

应助: 0 (幼儿园)
金币: 1456
帖子: 7
在线: 20.3小时
虫号: 255644
注册: 2006-05-28
性别: GG
专业: 原子和分子物理

[交流] 并行的一维fft fortran代码，速度还行

subroutine cfftserial ( n, x )
implicit none
integer (kind = 4) :: n, i
real (kind = 8)    :: pi, pid
complex (kind = 8) :: x(n, 2)

pi = acos ( -1.0d0 )
pid = -pi * 2.0d0
pid = pid / dble ( n )
!$omp parallel do default(shared) private(i)
do i = 0, n-1
   x(i+1, 1) = cmplx ( cos ( pid * dble ( i ) ), sin ( pid * dble ( i ) ), kind = 8 )
   x(i+1, 2) = conjg (  x(i+1, 1) )
end do
!$omp end parallel do

return
end subroutine cfftserial

subroutine ifftserial ( n, x )
implicit none
integer (kind = 4) :: n, i, j, k, m, it, is
integer (kind = 4) :: x(n)

k = 0
m = 1
do
   m = m * 2
   k = k + 1
   if ( n <= m ) exit
end do

!$omp parallel do default(shared) private(it, m, is, i, j)
do it = 0, n - 1
   m = it
   is = 0
   do i = 0, k - 1
      j = m / 2
      is = 2 * is + (m - 2 * j)
      m = j
   end do
   x(it + 1) = is + 1
end do
!$omp end parallel do

return
end subroutine ifftserial

subroutine purefour(n, isign, iserial, cserial, fft)
implicit none
integer (kind = 4) :: n, it, is, m, k, i, j, nv, l0, isign
integer (kind = 4) :: iserial(n)
real (kind = 8)    :: pi , pid, dn
complex (kind = 8) :: temp1, temp2, xsi
complex (kind = 8) :: fft(n), temp(n), cserial(n, 2)

dn = 1.0d0 / dble ( n )
k = 0
m = 1
do
   m = m * 2
   k = k + 1
   if ( n <= m ) exit
end do

!$omp parallel default(shared)

!$omp do private(i)
do i = 1, n
   temp(i) = fft(iserial(i))
end do
!$omp end do

!$omp do private(j, temp2)
do j = 0, n - 2, 2
   temp2 = temp(j + 1)
   temp(j + 1) = temp2 + temp(j + 2)
   temp(j + 2) = temp2 - temp(j + 2)
end do
!$omp end do

!$omp end parallel

if ( isign == 1) then
   is = 1
else
   is = 2
end if

m = n / 2
nv = 2

do l0 = k - 2, 0, -1
   m = m / 2
   nv = nv * 2
!$omp parallel do  &
!$omp default(private) &
!$omp shared(m, nv, temp, cserial, is)
   do it = 0, (m - 1) * nv, nv
      do j = 0, (nv / 2) - 1
         i = it + j + 1
         temp2 = cserial(m * j + 1, is) * temp(i + nv / 2)
         temp(i + nv / 2) = temp(i) - temp2
         temp(i) = temp(i) + temp2
      end do
   end do
!$omp end parallel do
end do

if ( isign == -1 ) then
!$omp parallel do &
!$omp default(shared)  private(i)
   do i = 1, n
      fft(i) = temp(i) * dn
   end do
!$omp end parallel do
else
!$omp parallel do &
!$omp default(shared)  private(i)
   do i = 1, n
      fft(i) = temp(i)
   end do
!$omp end parallel do
end if

return
end subroutine purefour

回复此楼