1)先上代码
!$omp parallel private(n, j, dx, dy, dz, r, a)
do n = 1, nsteps
!$omp do
do i = 0, nparticles - 1
x_tmp(i) = x(i) + vx(i) * dt
y_tmp(i) = y(i) + vy(i) * dt
z_tmp(i) = z(i) + vz(i) * dt
do j = 0, nparticles - 1
dx = x(j) - x(i)
dy = y(j) - y(i)
dz = z(j) - z(i)
R = sqrt(dx ** 2 + dy ** 2 + dz ** 2 + softening)
a = G * m(j) / R ** 3 * dt
vx(i) = vx(i) + a * dx
vy(i) = vy(i) + a * dy
vz(i) = vz(i) + a * dz
end do
end do
!$omp do
do i = 0, nparticles - 1
x(i) = x_tmp(i)
y(i) = y_tmp(i)
z(i) = z_tmp(i)
end do
!$omp master
print *, "n =", n
do i = 0, nparticles - 1
print *, i, ":", x(i), y(i), z(i)
end do
!$omp end master
end do
!$omp end parallel
2)运行结果
4 <-- [24, 16]
5 <-- [25, 17]
6 <-- [26, 18]
9 <-- [24]
10 <-- [25]
11 <-- [26]
13 <-- [9, 10, 11]
14 <-- [13]
16 <-- [14, 9]
17 <-- [14, 10]
18 <-- [14, 11]
24 <-- [4]
25 <-- [5]
26 <-- [6]
3)分析
看起来是对的。