版块导航: 正在加载中...

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

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

回复此楼

» 猜你喜欢

售SCI文章，我:8O.5.5.1O.54,科目全,可十急已经有11人回复
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急已经有9人回复
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急已经有11人回复
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急已经有12人回复
售SCI文章，我:8O.5.5.1O.54,科目全,可十急已经有9人回复
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急已经有9人回复
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急已经有11人回复
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急已经有11人回复
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急已经有8人回复
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急已经有8人回复

1楼 2011-12-26 17:00:35

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

找到一些相关的精华帖子，希望有用哦~

VASP并行计算时，提示“超过CPU时限”而退出计算。请问如何解决？已经有11人回复
Fortran 95／2003科学计算与工程对应的各章程序源代码已经有16人回复
经典书籍《数值算法大全（第3版）》：英文原版，含C、C++和Fortran源代码已经有303人回复
有没有人关注fortran语言的GPU并行呢？已经有12人回复
【求助】用最小二乘法求解线性方程组的Fortran代码已经有4人回复
【分享】C/C++常用算法设计方法,fortran其他语言编程也可以参考下已经有8人回复

点击这里搜索更多相关资源

科研从小木虫开始，人人为我，我为人人

相关版块跳转我要订阅楼主 panda0927 的主题更新