如何在每个时间步重新初始化 Fortran 函数

Question

我有一个程序如下：

Integer N,tmax
parameter (N=10,tmax=10)
real :: L,dt,noise,positions(tmax,N),positions0(N),s,force
integer :: param
force=0
L=100.0
dt=0.1
param=0
r=0
!initialization
do 1 i1=1,N
param=param+100
positions(1,i1)=(ran2(param)-0.5)*L
1 continue
positions0=positions(1,:)
!time dependent process
do 2 i2=1,tmax
    do 3 i3=1,N
    param=param+100
    force=0
    noise=0
    call routine_interaction(i3,N,positions0,force)
        noise=(ran2(param)-0.5)*L/10
        positions(i2,i3)=positions0(i3)+dt*(noise+force)
3   continue
positions0=positions(i2,:)
2 continue
print*, positions(1,:)
print*, positions(tmax,:)
end

子程序和函数：

      FUNCTION ran2(idum)
      INTEGER idum
!The ran2() subroutine from Numerical Recipes
!  (C) Copr. 1986-92 Numerical Recipes Software #$!5,5.){2p491&&k"15.  page 272
      END   


subroutine routine_interaction(ii,N,positions0,force)
integer N,ii
real :: r,s,force
real positions0(N)

        do 4 i4=1,N 
            if (ii.ne.i4) then
                r=abs(positions0(ii)-positions0(i4))
                s=sign(1.0,positions0(i4)-positions0(ii))
                force=force-(1/r**4-1/r**2)*s
            end if
4       continue
return
end

可以看到，在每个时间步都会调用子程序routine_interaction，这样在每个时间步，就有2次迭代和NxN=N²次计算。

fortran 中有没有一种方法可以在每个时间步定义一个函数，定义为：new_function(ii)=routine_interaction(ii,N,positions0,force) 并且会在 do 循环 3 中调用。这会导致在每个时间步进行 N+N 计算吗？

Answer 1

如果要使用这个算法，没有办法将力计算（routine_interaction）减少到O(N)时间。您在 positions0 中有 N 个元素，并且您正在计算每个元素与其他元素 (r=abs(positions(ii)-positions0(i4))) 之间的距离 r。这至少需要 N*(N-1)/2 次计算。

您可以使用三角循环稍微减少计算次数，例如作为

program p
  implicit none
    
  integer, parameter :: N=10, tmax=10
  real :: L,dt,noise,positions(tmax,N),positions0(N),force(N)
  integer :: param
  
  L=100.0
  dt=0.1
  param=0

  !initialization
  do i1=1,N
    param=param+100
    positions(1,i1)=(ran2(param)-0.5)*L
  enddo
  positions0=positions(1,:)
  
  !time dependent process
  do i2=1,tmax
    force = routine_interaction(N,positions0)
    do i3=1,N
      param=param+100
      noise=(ran2(param)-0.5)*L/10
      positions(i2,i3)=positions0(i3)+dt*(noise+force(i3))
    enddo
    positions0=positions(i2,:)
  enddo
  
  print*, positions(1,:)
  print*, positions(tmax,:)
end program

function routine_interaction(N,positions0) result (force)
  integer, intent(in) :: N
  real, intent(in) :: positions0(N)
  real :: r,s,force(N)
  integer :: i3,i4
  
  force = 0
  
  do i3=1,N
    do i4=1,i3-1
      r=abs(positions0(i3)-positions0(i4))
      s=sign(1.0,positions0(i4)-positions0(i3))
      force(i3)=force(i3)-(1/r**4-1/r**2)*s
      force(i4)=force(i4)+(1/r**4-1/r**2)*s
    enddo
  enddo
end function

如何在每个时间步重新初始化 Fortran 函数

How to reinitialize a Fortran function at each time step

performance

fortran

loops

function

gfortran