※ 引述《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) --

※ 发信站: 批踢踢实业坊(ptt.cc) ◆ From: 61.56.1.205