2次元配列matrixのループの順番

よく言われているように、FortranではC言語とは違い第1次元目が先にループするので、メインメモリからレジストやスタックに呼び出される際の効率の問題から(この辺、まだよくわかってません)、For文で計算する際は気をつけなければいけません。

それでどれくらい効率化するのかを調べるために、行列の掛け算を繰り返すサンプルプログラムを使用して時間を計測してみました。

プログラム内でやっていることは、要素が1000個の正方行列を作成し、10回掛け算を繰り返してans変数に格納しているだけです。出力すると、1000x1000の実数が出力されてコンソール画面上が埋め尽くされるので、コメントアウトしています。

記事の最後に、matrix.f90とmatrix_reverse.f90のソースコードを載せておきます。違いは、product_matrixサブルーチンなど、行列ループの回す順番です。よく見ると変数iとjが逆になっているはずです。

それぞれifortで最適化オプションは-O0にして計測しました。OSはScientific Linuxです。

matrix.f90 matrix_reverse.f90
real 2m14.443s 18m26.769s
user 2m14.237s 18m24.983s
sys 0m0.135s 0m0.988s

9倍くらい遅くなっています。


ちなみに、最適化オプションを-O3にした場合はこうなりました。

matrix.f90 matrix_reverse.f90
real 38.251s 40.885s
user 38.204s 40.842s
sys 0.022s 0.017s

コンパイル内でループの順番も最適化されているのでしょうか。もしくはメモリを節約したために1000x1000要素程度のmatrixは軽々と扱えるようになっているのでしょうか。


matrix.f90

program main
  implicit none

  integer(kind=8), parameter :: SEED_SELF = 1234567891234567
  integer(kind=8), parameter :: NELEM = 1000
  real(kind=8), allocatable, dimension(:,:) :: m1, ans
  real(kind=8) :: rd
  integer(kind=8) :: i
  integer(kind=8), parameter :: IMAX = 10

  allocate(m1(NELEM, NELEM))
  allocate(ans(NELEM, NELEM))

  ! zero clear
  write(*,*) "## zero clear"
  call init_matrix(m1)
  call init_matrix(ans)

  ! set random number
  write(*,*) "## set rundom number"
  call setrand_matrix(m1, SEED_SELF)
  ans(:,:) = m1(:,:)

  do i = 1, IMAX
    ! calculate
    write(*,*) "## calculate product"
    call product_matrix(ans, m1, ans)
    !call print_matrix(ans)
  enddo

  ! output
  write(*,*) "## output matrix : ",IMAX
  !call print_matrix(ans)

  deallocate(m1)
  deallocate(ans)

  contains

  subroutine init_matrix(mat)
    implicit none
  
    real(kind=8), dimension(:,:), intent(inout) :: mat
  
    integer(kind=8) :: i, j
    integer(kind=8) :: nelem
  
    nelem = size(mat, 1)
  
    do i = 1, nelem
      do j = 1, nelem
        mat(j, i) = 0.0d0
      enddo
    enddo
    
  
  end subroutine init_matrix
  
  subroutine setrand_matrix(mat, iseed)
    implicit none
  
    real(kind=8), dimension(:,:), intent(inout) :: mat
    integer(kind=8), intent(in) :: iseed
  
    integer(kind=8), allocatable, dimension(:) :: seed
    integer(kind=8) :: nrand
    integer(kind=8) :: i, j
    integer(kind=8) :: nelem
  
    seed = iseed
    call random_seed(put=seed)
  
    nelem = size(mat, 1)
  
    do i = 1, nelem
      do j = 1, nelem
        call random_number(mat(j, i))
      enddo
    enddo
    
  
  end subroutine setrand_matrix
  
  subroutine print_matrix(mat)
    implicit none
  
    real(kind=8), dimension(:,:), intent(in) :: mat
  
    integer(kind=8) :: i, j
    integer(kind=8) :: nelem
  
    nelem = size(mat, 1)
  
    write(*,fmt='(8x)', advance='no')
    do j = 1, nelem
      write(*,fmt='(1x, i12)', advance='no') j
    enddo
    write(*,*)
  
    do i = 1, nelem
      write(*,fmt='(1x, i7)', advance='no') i
      do j = 1, nelem
        write(*,fmt='(1x, e12.5)', advance='no') mat(i, j)
      enddo
      write(*,*)
    enddo
    
  
  end subroutine print_matrix
  
  subroutine product_matrix(m1, m2, ans)
    implicit none
  
    real(kind=8), dimension(:,:), intent(in) :: m1, m2
    real(kind=8), dimension(:,:), intent(out) :: ans
  
    integer(kind=8) :: i, j, k
    integer(kind=8) :: nelem
  
    nelem = size(m1, 1)
  
    do i = 1, nelem
      do j = 1, nelem
        do k = 1, nelem
          ans(j, i) = m1(j, k) * m2(k, i)
        enddo
      enddo
    enddo
  
  end subroutine product_matrix
end program main

matrix_reverse.f90

program main
  implicit none

  integer(kind=8), parameter :: SEED_SELF = 1234567891234567
  integer(kind=8), parameter :: NELEM = 2000
  real(kind=8), allocatable, dimension(:,:) :: m1, ans
  real(kind=8) :: rd
  integer(kind=8) :: i
  integer(kind=8), parameter :: IMAX = 10

  allocate(m1(NELEM, NELEM))
  allocate(ans(NELEM, NELEM))

  ! zero clear
  write(*,*) "## zero clear"
  call init_matrix(m1)
  call init_matrix(ans)

  ! set random number
  write(*,*) "## set rundom number"
  call setrand_matrix(m1, SEED_SELF)
  ans(:,:) = m1(:,:)

  do i = 1, IMAX
    ! calculate
    write(*,*) "## calculate product"
    call product_matrix(ans, m1, ans)
    !call print_matrix(ans)
  enddo

  ! output
  write(*,*) "## output matrix : ",IMAX
  !call print_matrix(ans)

  deallocate(m1)
  deallocate(ans)

  contains

  subroutine init_matrix(mat)
    implicit none
  
    real(kind=8), dimension(:,:), intent(inout) :: mat
  
    integer(kind=8) :: i, j
    integer(kind=8) :: nelem
  
    nelem = size(mat, 1)
  
    do i = 1, nelem
      do j = 1, nelem
        mat(i, j) = 0.0d0
      enddo
    enddo
    
  
  end subroutine init_matrix
  
  subroutine setrand_matrix(mat, iseed)
    implicit none
  
    real(kind=8), dimension(:,:), intent(inout) :: mat
    integer(kind=8), intent(in) :: iseed
  
    integer(kind=8), allocatable, dimension(:) :: seed
    integer(kind=8) :: nrand
    integer(kind=8) :: i, j
    integer(kind=8) :: nelem
  
    seed = iseed
    call random_seed(put=seed)
  
    nelem = size(mat, 1)
  
    do i = 1, nelem
      do j = 1, nelem
        call random_number(mat(i, j))
      enddo
    enddo
    
  
  end subroutine setrand_matrix
  
  subroutine print_matrix(mat)
    implicit none
  
    real(kind=8), dimension(:,:), intent(in) :: mat
  
    integer(kind=8) :: i, j
    integer(kind=8) :: nelem
  
    nelem = size(mat, 1)
  
    write(*,fmt='(8x)', advance='no')
    do j = 1, nelem
      write(*,fmt='(1x, i12)', advance='no') j
    enddo
    write(*,*)
  
    do i = 1, nelem
      write(*,fmt='(1x, i7)', advance='no') i
      do j = 1, nelem
        write(*,fmt='(1x, e12.5)', advance='no') mat(i, j)
      enddo
      write(*,*)
    enddo
    
  
  end subroutine print_matrix
  
  subroutine product_matrix(m1, m2, ans)
    implicit none
  
    real(kind=8), dimension(:,:), intent(in) :: m1, m2
    real(kind=8), dimension(:,:), intent(out) :: ans
  
    integer(kind=8) :: i, j, k
    integer(kind=8) :: nelem
  
    nelem = size(m1, 1)
  
    do i = 1, nelem
      do j = 1, nelem
        do k = 1, nelem
          ans(i, j) = ans(i, j) + m1(i, k) * m2(k, j)
        enddo
      enddo
    enddo
  
  end subroutine product_matrix
end program main