※ 引述《DMFC (sole)》之铭言：
: ※ 引述《Yagyu (近在眼前)》之铭言：
: : do i=1,nx ; do j=1,ny ; do k=1,nz
: : csum=sum(coeff(1:np)*cdexp)/Vol
: : csumout(i,j,k)=csum
: : end do ; end do ; end do
: : sum(coeff(1:np)*cdexp)
: : 请问这边是要将整个 coeff(1:23000) 乘上 cdexp 再做 sum 吗?
: : 会这么问是因为我不确定你未贴上程式码部分是否还有 np 的 loop
: : 如果没有 np loop, 同时 cdexp 只是个常数,不会随 nx,ny,nz 变动
: : 那可以试着将这部份移出循环外, 这样省得每次 loop 都要重算一遍
: : 如果有 np loop, 那请无视我的想法
: : 另外输出的部份
: : do i=1,nx
: : do j=1,ny
: : do k=1,nz
: : csumout(i,j,k)=sum(coeff(1:np)*cdexp)/Vol
: : end do
: : end do
: : end do
: : do k=1,nz
: : do j=1,ny
: : do i=1,nx
: : csumout(i,j,k)=sum(coeff(1:np)*cdexp)/Vol
: : end do
: : end do
: : end do
: : 两者差异 请参考彭国伦先生的fortran工具书 应该是在阵列章节中的多维阵列那边
: : 没记错的话 是跟内存存放资料方式有关 这边变动我想绝对有帮助
: 谢谢
: 没错~顺序影响很大
: 不过这是因为我手误
: 由于我没法直接贴SOURCE CODE
: 所以是用手KEYIN
: LOOP顺序我贴错了
: 且~很无奈的
: 我那个 cdexp 是与 i,j,k 有关
: cdexp 非 常数~无法提出
: 我再贴一次完整的CODE
: do iz=0,ngrid(3)-1 ; do iy=0,ngrid(2)-1 ; do ix=0,ngrid(1)-1
: xyz(1) = dble(ix)/dble(ngrid(1))
: xyz(2) = dble(iy)/dble(ngrid(2))
: xyz(3) = dble(iz)/dble(ngrid(3))
: atmp = pi2 * (wkiG1*xyz(1) + wkiG2*xyz(2) + wkiG3*xyz(3))
: csumout(ix,iy,iz) = sum(coeff(1:nplane)*cdexp(atmp(1:nplane)))/dsqrt(Vol)
^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^
1D array 1D array
^^^^^^^^^^^^^^^^^^^^^
atmp里面的element全部作cdexp
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
再作array相乘
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
再sum
: xyz(1) = dble(ix)/dble(ngrid(1))
: xyz(2) = dble(iy)/dble(ngrid(2))
: xyz(3) = dble(iz)/dble(ngrid(3))
这边浪费很多加减乘除的处理时间...推文中讲过...不再累述....
程式码无全数贴出故无法得知中间少了啥....囧
理论上atmp & coeff这两个array要先填完...接下去作才有意义...
建议你在atmp = ......这一行用Euler formular拆解re part and im part
并先算完cdexp后...
再到下一行直接用dot函数处理并填到csumout里头....
最后3层loop结束后再加一行 csumout=csumout/dsqrt(Vol)
: end do ; end do ; end do
: ngrid(1) ngrid(2) ngrid(3) 都是常数
: atmp, wkiG1, wkiG2, wkiG3 都是维度23000(nplane)的大矩阵
: Vol, pi2 是常数
: cdexp 是 fortran 内有默认的function
: 意思是对 double precision 的 complex 取 exp
: (atmp 是double precision 的 complex)