FORTRAN:通过指针矩阵访问数组,性能良好



我在这里遇到了使用指针的问题。在我这么做之前,我有一个性能问题。假设有一个像这样的二维矩阵:

0.0  0.0  0.0.....
0.0  0.7  0.5.....
0.0  0.5  0.8.....
0.0  0.3  0.8.....

我需要计算这个东西的梯度。因此,对于每个数字,我需要这个数字以及这个2D矩阵的所有4个最近邻居。此外,第一行和最后一行及列为0。

现在我有两种方法:

  1. 直接制作这样一个NxN矩阵并计算梯度。完全按照描述。这里的内存使用是NxNxreal*8,循环从计算(2,2)元素开始,然后计算(2,3)。。。

  2. 制作一个(N-2)x(N-1)+1数组和一个NxN指针矩阵(此时使用类型)。数组的(N-2)x(N-2个)元素将存储除边界上的0.0s之外的数字。的最后一个元素是0.0。对于指针矩阵,边界上的所有元素都将指向数组的最后一个元素0.0。其他指针应该指向它们应该指向的地方。

性能问题来了,因为我处理的矩阵可能非常巨大,也可能是3D的。

对于方法1,没有什么好说的,因为它只是一个直接的方法。

对于方法2,我想知道编译器是否能够正确处理这个问题。根据我目前的理解,每个FORTRAN指针都像一个结构。如果是这种情况,FORTRAN指针比c指针慢,因为它不仅仅是一个简单的去引用我还想知道指针的类型扭曲是否会降低性能(制作指针矩阵需要这种扭曲)。有没有特别的原因让我放弃方法2,因为它应该更慢?

让我们以windows上的IVF、Linux上的gfortran和ifort为例。因为它可以依赖于编译器。

更新:感谢Stefan的代码。我自己也写了。

program stencil
    implicit none
    type pp
        real*8, pointer :: ptr
    endtype pp
    type(pp), allocatable :: parray(:,:)
    real*8, allocatable, target :: array(:)
    real*8, allocatable :: grad(:,:,:), direct(:,:)
    integer, parameter :: n = 5000
    integer :: i, j
    integer :: clock_rate, clock_start, clock_stop
    allocate(array(n**2+1))
    allocate(parray(0:n+1, 0:n+1))
    allocate(grad(2, n, n))
    call random_number(array)
    array(n**2+1) = 0
    do i = 0, n + 1
        parray(0,i)%ptr => array(n**2+1)
        parray(n+1,i)%ptr => array(n**2+1)
        parray(i,0)%ptr => array(n**2+1)
        parray(i,n+1)%ptr => array(n**2+1)
    enddo
    do i = 1, n
        do j = 1, n
            parray(i,j)%ptr => array((i-1) * n + j)
        enddo
    enddo
    !now stencil
    call system_clock(count_rate=clock_rate)
    call system_clock(count=clock_start)
    do j = 1, n
        do i = 1, n
            grad(1, i, j) = (parray(i + 1,j)%ptr - parray(i - 1,j)%ptr)/2.D0
            grad(2, i, j) = (parray(i,j + 1)%ptr - parray(i,j - 1)%ptr)/2.D0
        enddo
    enddo
    call system_clock(count=clock_stop)
    print *, "pointer, time cost= ", real(clock_stop-clock_start)/real(clock_rate)
    deallocate(array)
    deallocate(parray)
    allocate(direct(0:n+1, 0:n+1))
    call random_number(direct)
    do i = 0, n + 1
        direct(0,i) = 0
        direct(n+1,i) = 0
        direct(i,0) = 0
        direct(i,n+1) = 0
    enddo
    !now stencil directly
    call system_clock(count_rate=clock_rate)
    call system_clock(count=clock_start)
    do j = 1, n
        do i = 1, n
            grad(1, i, j) = (direct(i + 1,j) - direct(i - 1,j))/2.D0
            grad(2, i, j) = (direct(i,j + 1) - direct(i,j - 1))/2.D0
        enddo
    enddo
    call system_clock(count=clock_stop)
    print *, "direct, time cost= ", real(clock_stop-clock_start)/real(clock_rate)
endprogram stencil

结果(o0):

指针,时间成本=2.170000

直接,时间成本=1.127000

结果(o2):

指针,时间成本=0.5110000

直接,时间成本=9.4999999E-02

所以FORTRAN指针要慢得多。Stefan早些时候已经指出了这一点。现在我想知道是否还有改进的余地。据我所知,如果我用c来做,应该不会有太大的区别。

首先,我必须道歉,因为我误解了指针在Fortran中的工作方式。。。


最后,我对这个话题很感兴趣,于是我自己做了一个测试。它基于一个数组,该数组的周围有一个零。

声明:

real, dimension(:,:), allocatable, target :: array
real, dimension(:,:,:), allocatable :: res
real, dimension(:,:), pointer :: p1, p2, p3, p4
allocate(array(0:n+1, 0:n+1), source=0.)
allocate(res(n,n,2), source=0.)

现在的方法:

循环:

do j = 1, n
    do i = 1, n
        res(i,j,1) = array(i+1,j) - array(i-1,j)
        res(i,j,2) = array(i,j+1) - array(i,j-1)
    end do
end do

阵列分配:

res(:,:,1) = array(2:n+1,1:n) - array(0:n-1,1:n)
res(:,:,2) = array(1:n,2:n+1) - array(1:n,0:n-1)

指针:

p1 => array(0:n-1,1:n)
p2 => array(1:n,2:n+1)
p3 => array(2:n+1,1:n)
p4 => array(1:n,0:n-1)
res(:,:,1) = p3 - p1
res(:,:,2) = p2 - p4

虽然最后两种方法确实依赖于额外的零层,但循环可以引入一些条件来处理这些问题。

时间安排很有趣:

loops:     0.17528710301849060
array:     0.21127231500577182
pointers:  0.21367537401965819

虽然数组和指针分配产生的时间大致相同,但循环构造(注意循环顺序!这是5的倍数!!)是最快的方法。


更新:我试图从你的代码中挤出一点性能,发现了一件小事。您的代码在0.95s0.30s(使用n = 10000)中使用-O2执行。

转换矩阵以获得更线性的内存访问,指针部分的运行时间为0.50s

parray(i,j)%ptr => array((j-1) * n + i)

IMHO,问题是缺少关于指针的信息,这禁止额外的优化。使用-O3 -fopt-info-missed,您会收到关于未知对齐和非连续访问的投诉。与我的结果相比,额外的因素2应该源于这样一个事实,即您使用的是双精度,而我的代码是以单精度编写的。。。

我认为Stefan的答案是最好的答案。但我个人想为讨论和我自己的问题做一个结论。

  1. 根据Vladimir的说法,FORTRAN指针不同于C指针。FORTRAN标准似乎旨在使数组指针成为数组的"子对象"。因此,与C中的情况不同,FORTRAN中的"指针数组"几乎毫无意义。有关使用FORTRAN指针的详细信息,请阅读Stefan的代码。此外,FORTRAN中的"指针数组"是可能的,但它的低性能绝对不是一个简单的解引用。

  2. 使用循环展开的直接访问可以提高计算性能。请在Stefan的代码中找到详细信息。当使用直接访问时,根据Stefan的评论,编译器优化可以做得更好。我认为这就是为什么人们不使用指针来解决模具问题的原因。

  3. 使用指针处理模具的想法是降低内存成本,并使边界条件非常灵活。但它暂时不是性能的选择。主要原因是在不了解指针模式的情况下,无法执行"非连续"内存访问和编译器优化。

关于FORTRAN指针,请参阅Stefan的回答。

最新更新