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