[問題] 迴圈拆開計算後與原本值有微小誤差

看板Fortran作者boa85391 ( )時間9年前 (2015/08/12 21:27)推噓2(2推 0噓 6→)

留言8則, 2人參與討論串1/2 (看更多)

do i=1, 3999, 1 do j=i+1, 4000, 1 rx=bead(j)%x_p - bead(i)%x_p ry=bead(j)%y_p - bead(i)%y_p rz=bead(j)%z_p - bead(i)%z_p rpq= dsqrt(rx**2+ry**2+rz**2) if (rpq<=r_c) then if (bead(i)%bead_number==1 .and. bead(j)%bead_number==1) then rep_term=dabs(24.d0*eps_hh*(2.d0*sigma**12*rpq**(-14)& -sigma**6*rpq**(-8))) bead(i)%x_frep=bead(i)%x_frep - rep_term * rx bead(i)%y_frep=bead(i)%y_frep - rep_term * ry bead(i)%z_frep=bead(i)%z_frep - rep_term * rz bead(j)%x_frep=bead(j)%x_frep + rep_term * rx bead(j)%y_frep=bead(j)%y_frep + rep_term * ry bead(j)%z_frep=bead(j)%z_frep + rep_term * rz else if (bead(i)%bead_number/=1 .and. bead(j)%bead_number/=1) then rep_term=dabs(24.d0*eps_tt*(2.d0*sigma**12*rpq**(-14)& -sigma**6*rpq**(-8))) bead(i)%x_frep=bead(i)%x_frep - rep_term * rx bead(i)%y_frep=bead(i)%y_frep - rep_term * ry bead(i)%z_frep=bead(i)%z_frep - rep_term * rz bead(j)%x_frep=bead(j)%x_frep + rep_term * rx bead(j)%y_frep=bead(j)%y_frep + rep_term * ry bead(j)%z_frep=bead(j)%z_frep + rep_term * rz if (bead(i)%Lipid_number/=bead(j)%Lipid_number) then bead(i)%x_fatt=bead(i)%x_fatt - rep_term * rx bead(i)%y_fatt=bead(i)%y_fatt - rep_term * ry bead(i)%z_fatt=bead(i)%z_fatt - rep_term * rz bead(j)%x_fatt=bead(j)%x_fatt + rep_term * rx bead(j)%y_fatt=bead(j)%y_fatt + rep_term * ry bead(j)%z_fatt=bead(j)%z_fatt + rep_term * rz end if else rep_term=dabs(24.d0*eps_ht*(2.d0*sigma**12*rpq**(-14)& -sigma**6*rpq**(-8))) bead(i)%x_frep=bead(i)%x_frep - rep_term * rx bead(i)%y_frep=bead(i)%y_frep - rep_term * ry bead(i)%z_frep=bead(i)%z_frep - rep_term * rz bead(j)%x_frep=bead(j)%x_frep + rep_term * rx bead(j)%y_frep=bead(j)%y_frep + rep_term * ry bead(j)%z_frep=bead(j)%z_frep + rep_term * rz end if else if (r_c<=rpq .and. rpq<=(r_c+w_tt) .and. & bead(i)%Lipid_number/=bead(j)%Lipid_number .and. & bead(i)%bead_number/=1 .and. bead(j)%bead_number/=1) then att_term=dabs(eps_tt*pi/w_tt/rpq*dcos(pi*(rpq-r_c)/2.d0/w_tt)* & dsin(pi*(rpq-r_c)/2.d0/w_tt)) bead(i)%x_fatt=bead(i)%x_fatt + att_term * rx bead(i)%y_fatt=bead(i)%y_fatt + att_term * ry bead(i)%z_fatt=bead(i)%z_fatt + att_term * rz bead(j)%x_fatt=bead(j)%x_fatt - att_term * rx bead(j)%y_fatt=bead(j)%y_fatt - att_term * ry bead(j)%z_fatt=bead(j)%z_fatt - att_term * rz end if end do end do 主要是有4000顆珠子，計算彼此間距離一條件算作用力，rep是排斥力，att是吸引力，除了if判定式裡面的number存成integer，其他都是real(8) 為了方便使用openMP，所以我把它寫成兩個與上式一樣的雙層迴圈，一個只算i的部分，一個只算j的部分，可是最後得到的值有些會在小數點後十幾位有誤差嘗試輸出每一步的rep_term去比較，存到小數點後15位都還是沒什麼問題想請問一下這是合理的誤差產生嗎?還是我的寫法會結構上有問題，已經在這部分弄了快兩天，因為是直觀的同樣迴圈算兩次，還是無法去忽略這些誤差。 -- ※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.112.23.59 ※ 文章網址: https://www.ptt.cc/bbs/Fortran/M.1439386065.A.BAC.html