よく言われているように、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