fortran - 为什么按元素矩阵行交换比 Fortran 中按数组行交换更有效？

Question

我有一些执行矩阵行交换的 Fortran 代码。在遗留代码中，它写成

do J = 1,N
   save_val  = B(IROW,J)
   B(IROW,J) = B(JCOL,J)
   B(JCOL,J) = save_val
end do

这会将 IROW 行与 JCOL 行进行交换(IROW 和 JCOL 是整数)。但是，该代码块的功能并不直观。在我看来，这样写会更直观，或者至少有助于提高可读性:

save_row  = B(IROW,:)
B(IROW,:) = B(JCOL,:)
B(JCOL,:) = save_row

(更清楚的是行正在移动)。

从附图中可以清楚地看出，循环方法相对于数组操作提供了更好的性能。为什么是这样？是因为当数组中的元素数量变大时，这会变成一个内存受限的过程吗？ (即数组会被“分块”)还是其他原因？

编译为 gfortran -O3 test.f95。添加标志 fstack-arrays 没有产生显着差异。

program test

  implicit none

  integer  :: N
  integer  :: M
  integer  :: loop_max = 1e7
  integer  :: i                    ! loop index
  real     :: t1, t2
  real     :: time_loop, time_array, time_sub_loop, time_sub_array

  real, dimension(:, :), allocatable   :: B
  real, dimension(:)   , allocatable   :: save_row

  real :: save_val
  integer :: IROW, J, JCOL

  character(*), parameter :: format_header = '(A5, 1X, 4(A12,1X))'
  character(*), parameter :: format_data = '(I5, 1X, 4(ES12.5, 1X))'


  open(1, file = 'TimingRowExchange.txt', status = 'unknown')
  write(1, format_header) 'N', 't_loop', 't_array', 't_sub_loop', 't_sub_array'

  do N = 1, 100
    M = N + 1
    allocate(B(N,N), save_row(M))
    call random_number(B)

    JCOL = 1
    IROW = 3

    call CPU_time(t1)
    do i = 1, loop_max
      do J = 1,N
        save_val  = B(IROW,J)
        B(IROW,J) = B(JCOL,J)
        B(JCOL,J) = save_val
      end do
    end do
    call CPU_time(t2)
    time_loop = t2 - t1
    ! write ( *, * ) 'Using Loop =', t2 - t1


    call CPU_time(t1)
    do i = 1, loop_max
        save_row(1:N) = B(IROW,:)
        B(IROW,:)     = B(JCOL,:)
        B(JCOL,:)     = save_row(1:N)
    end do
    call CPU_time(t2)
    time_array = t2 - t1
    ! write ( *, * ) 'Using Array =', t2 - t1

    call CPU_time(t1)
    do i = 1, loop_max
      call exchange_rows_loop(B, JCOL, IROW)
    end do
    call CPU_time(t2)
    time_sub_loop = t2 - t1
    ! write ( *, * ) 'Loop Subrout =', t2 - t1


    call CPU_time(t1)
    do i = 1, loop_max
      call exchange_rows_array(B, JCOL, IROW)
    end do
    call CPU_time(t2)
    time_sub_array = t2 - t1
    ! write ( *, * ) 'Array Subrout =', t2 - t1

    deallocate(B, save_row)
    write(1, format_data) N, time_loop, time_array, time_sub_loop, time_sub_array
  end do


contains


  subroutine print_mat(A)
    implicit none
    real, dimension(:,:), intent(in) :: A
    integer :: n

    n = size(A,1) ! # of rows

    do i = 1,n
      print*, A(i,:)
    end do
    print*,

  end subroutine print_mat



  subroutine exchange_rows_loop(A, row1, row2)
    implicit none
    real, dimension(:,:), intent(in out) :: A
    integer,              intent(in)     :: row1, row2

    integer :: J
    real :: save_val

    do J = 1, size(A,1)
      save_val  = A(row1,J)
      A(row1,J) = A(row2,J)
      A(row2,J) = save_val
    end do

  end subroutine exchange_rows_loop



  subroutine exchange_rows_array(A, row1, row2)
    implicit none
    real, dimension(:,:), intent(in out) :: A
    integer,              intent(in)     :: row1, row2

    real, dimension(size(A,1))           :: save_row

    save_row  = A(row1,:)
    A(row1,:) = A(row2,:)
    A(row2,:) = save_row

  end subroutine exchange_rows_array


end program test

Answer 1

我对 Fortran 哲学(优势)的理解是，该语言应该帮助用户专注于科学，同时处理大多数与计算机相关的事情，例如速度优化、垃圾收集等。

通过 pure/elemental 函数和子例程的函数式编程风格是恕我直言，已被引入但未得到充分利用的最伟大的工具之一，因为它使代码更清晰、更简单且更丰富健壮。

所以我又添加了一个带有 elemental 交换例程的测试:

  subroutine exchange_rows_elemental(A, row1, row2)
    implicit none
    real, dimension(:,:), intent(in out) :: A
    integer,              intent(in)     :: row1, row2
    call swap(A(row1,:),A(row2,:))
  end subroutine exchange_rows_elemental  

  elemental subroutine swap(a,b)
     real, intent(inout) :: a,b
     real :: save_val
     save_val = a
     a = b
     b = save_val
  end subroutine swap  

主要是:

call CPU_time(t1)
do i = 1, loop_max
  call exchange_rows_elemental(B, JCOL, IROW)
end do
call CPU_time(t2)
time_elemental = t2 - t1
! write ( *, * ) 'Elemental =', t2 - t1

这是我在 Windows 上使用 gfortran 9.2.0 得到的结果:

elemental 版本几乎与最快的循环版本一样快，但它可能以矢量化方式运行。我确定在这种情况下，编译器可能正在内联 swap 例程(如果它在另一个文件中，它可能无法这样做)，但仍然告诉编译器 swap 例程可以矢量化可能有助于它实现最佳性能。我喜欢它，因为它是一种很好的方式，可以充分利用编译器优化，而不会用嵌套循环和循环变量使源代码困惑。