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