使用广播的概念获取 Fortran 数组中 2D 空间中点的距离 (Python)

Question

我是 Fortran 新手，我已经很难理解 Python 中数组广播的概念，所以在 Fortran 中实现它对我来说更加困难

Fortran 代码：

program test
    implicit none
    integer,parameter::N=6,t_size=500
    real,dimension(t_size,N,2)::array1
    
    contains

    pure function f(a) result(dr)
    real::dr(N,N,2)
    real,intent(in)::a(N,2)
    real::b(N,N,2)
    real::c(N,N,2)
    b=spread(a,1,N)
    c=spread(a,2,N)
    dr=c-b

    end function f

end program

现在 N 是点数，t_size 只是不同时间步长的数量。我想出了这个函数，它在两个不同的维度中使用 spread 来创建 NxNx2 数组。

现在我想到使用像

r=f(array1(1,:,:)

这样的线来获得一个数组，该数组保存每2个点的空间坐标的所有差异。

我已经用 Python 编写了执行此操作的代码（取自 Python 物理教科书）

r = np.empty((t.size, N, 2))
r[0] = r0

def f(r):
    dr=r.reshape(N,1,2)-r

我稍后可以写的地方，例如

f(r[i]

。 （在本例中，我留下了 r[0] = r0 行，因为它表明给出了初始条件 - 稍后我计划在 Fortran 中使用 random_number 子例程来执行此操作。）

现在我希望我的问题清楚了。如果有人有更好的想法（我确信有）在 Fortran 中实现广播，请告诉我。

请对 Fortran 新手和一般编程新手有一点耐心 - 提前感谢您的回复。

我已经尝试过使用 random_number 子例程并且它有效，但我无法检查输出是否为真。

Answer 1

现在我有更多时间，我会将我的评论变成答案。

Fortran 是一种编译语言。 Python 是一种解释型语言。这导致了许多差异，但最大的差异之一是 Fortran 循环速度很快，因为它们已转换为优化的机器代码，并且通常是表达您想要执行的操作的首选方式。然而，在 Python 中，循环很慢，因为它们必须（重新）解释每一步 - 因此鼓励使用（在我看来）神秘的单行程序从 python 库调用优化函数。 Fortran 中的简单操作很少需要这种“向量化”——（在我看来）只需用循环的简单方式编写它即可。

当我在这里时，我还会提到，如果您确实想要性能，那么您的数组索引的方式是错误的 - 您希望最常更改的索引（上面的 1:2 ）是 first 索引。这最终与 Fortran 在内存中组织数组的方式有关，如果有兴趣，要研究的关键术语是“列主”

最后我要提一下，Fortran 的方式通常是使用

Subroutine

而不是

Function

。这主要是风格上的，但确实导致了一些简化，我不会在这里讨论。

无论如何，将以上内容放在一起，我会写成如下所示

ijb@ijb-Latitude-5410:~/work/stack$ cat not_python.f90
Program test

  Implicit None

  Real, Dimension( :, :, : ), Allocatable :: array1
  Real, Dimension( :, :, : ), Allocatable :: dr
  
  Integer :: d, N, t_size
  Integer :: t_step
  Integer :: id, iN, it
  
  Write( *, * ) 'N, t_size ?'
  Read ( *, * )  N, t_size

  d = 2

  Allocate( array1( 1:d, 1:N, 1:t_size ) )

  Do it = 1, Size( array1, Dim = 3 )
     Do iN = 1, Size( array1, Dim = 2 )
        Do id = 1, Size( array1, Dim = 1 )
           array1( id, iN, it ) = id + 10.0 * in + 100.0 * it
        End Do
     End Do
  End Do

  Allocate( dr( 1:d, 1:n, 1:n ) )

  t_step = 1
  Call f( array1( :, :, t_step ), dr )

  Call printit( dr )
  

Contains

  Pure Subroutine f( a, dr )

    Implicit None

    Real, Dimension(    :, : ), Intent( In    ) :: a
    Real, Dimension( :, :, : ), Intent(   Out ) :: dr

    Integer :: N
    Integer :: i, j
    
    N = Size( a, Dim = 2 )

    ! Could take advantage of symmetry here, dr( j, i ) = - dr( i, j )
    ! but probably gaines you little 
    Do j = 1, N
       Do i = 1, N
          dr( :, i, j ) = a( :, i ) - a( :, j )
       End Do
    End Do

  End Subroutine f

  Subroutine printit( dr )

    Use iso_fortran_env, Only : stdout => output_unit

    Implicit None
    
    Real, Dimension( :, :, : ), Intent( In ) :: dr

    Integer :: it, iN, id
    
    Do it = 1, Size( dr, Dim = 3 )
       Do iN = 1, Size( dr, Dim = 2 )
          Do id = 1, Size( dr, Dim = 1 )
             Write( stdout, '( 3( i2, 1x ), 6x, f5.1 )' ) &
                  id, iN, it, dr( id, iN, it )
          End Do
       End Do
    End Do

  End Subroutine printit

End Program test
ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -O -std=f2018 not_python.f90 -o not_python
ijb@ijb-Latitude-5410:~/work/stack$ ./not_python 
 N, t_size ?
3 2
 1  1  1         0.0
 2  1  1         0.0
 1  2  1        10.0
 2  2  1        10.0
 1  3  1        20.0
 2  3  1        20.0
 1  1  2       -10.0
 2  1  2       -10.0
 1  2  2         0.0
 2  2  2         0.0
 1  3  2        10.0
 2  3  2        10.0
 1  1  3       -20.0
 2  1  3       -20.0
 1  2  3       -10.0
 2  2  3       -10.0
 1  3  3         0.0
 2  3  3         0.0