2次元配列: Sample 2: マトリックス・ベクトルの掛け算(FORTRAN)

高速化プログラミング   
トップ  >  メモリジャンプと高速化  >  2次元配列  >  Sample 2: マトリックス・ベクトルの掛け算(FORTRAN)

2次元配列: Sample 2: マトリックス・ベクトルの掛け算(FORTRAN)

言語の変更:   C版   JScript版

■ 概要

このサンプルではマトリックスとベクトルの掛け算を行います。

c=Ab           ... 式(1)

通常、マトリックスAijを2次元配列A[i][j]としてコーディングしますが、場合によっては列と行を入れ替えてA[j][i]にすることもできます。ベクトルcさえ求まればどちらも採用してもいいですが、問題は計算時間です。実際にプログラムを見て、実行してみましょう。

下のCode 1ではA[j][i]、Code 2ではA[i][j]としてコーディングしてマトリックスとベクトルの掛け算をします。マトリックス(ベクトル)のサイズnは1000から30000までの30個の値を用います。各nに対して計算にかかった時間を測定し出力します。

■ ソースコード


  ◆ Code 1   ◆ Code 2
 
c Program to measure time to carry out matrix-vector product.
c c(i) = sum(a(i,j)*b(j),j=1,...,n) for i=1,...,n
      program main
      implicit none
      integer i, j, k, n
      integer, allocatable :: a(:,:), b(:), c(:)
      integer*8 time0, time1, dtime
      real*8 time
c
      write(*,*)"Matrix size    Elapsed time [sec]"
      do k=1,30
c Matrix size
        n = k * 1000

c Allocation
        allocate(a(n,n), b(n), c(n))

c Initialization
        do i=1,n
          do j=1,n
            a(j,i) = 1
          end do
          b(i) = 1
        end do

c Start time
        call system_clock(time0)

c Main calculation: matrix-vector product
        do i=1,n
          c(i) = 0
          do j=1,n
            c(i) = c(i) + a(i,j) * b(j)
          end do
        end do

c Finish time
        call system_clock(time1, dtime)

c Output time
        time = 1d0*(time1-time0)/dtime
        write(*,"(i12,f16.7)")n, time

c Deallocation
        deallocate(a, b, c)
      end do

      end program

    
 
c Program to measure time to carry out matrix-vector product.
c c(i) = sum(a(j,i)*b(j),j=1,...,n) for i=1,...,n
      program main
      implicit none
      integer i, j, k, n
      integer, allocatable :: a(:,:), b(:), c(:)
      integer*8 time0, time1, dtime
      real*8 time
c
      write(*,*)"Matrix size    Elapsed time [sec]"
      do k=1,30
c Matrix size
        n = k * 1000

c Allocation
        allocate(a(n,n), b(n), c(n))

c Initialization
        do i=1,n
          do j=1,n
            a(j,i) = 1
          end do
          b(i) = 1
        end do

c Start time
        call system_clock(time0)

c Main calculation: matrix-vector product
        do i=1,n
          c(i) = 0
          do j=1,n
            c(i) = c(i) + a(j,i) * b(j)
          end do
        end do

c Finish time
        call system_clock(time1, dtime)

c Output time
        time = 1d0*(time1-time0)/dtime
        write(*,"(i12,f16.7)")n, time

c Deallocation
        deallocate(a, b, c)
      end do

      end program

    

■ 計算時間の測定結果

Code 1Code 2を実行したときの計算時間をそれぞれ青線と赤線で図1に示します。そして両者の比を緑線で示します。

測定時間
図1 測定時間

■ 考察

FORTRAN言語では2次元配列を列単位で1次元配列に組みなおされるので、内側のループが行に関連付けられたCode 1は図2のように飛び飛びしたメモリアクセスになります。一方、内側のループが行に関連付けられたCode 2は図3のように連続したメモリアクセスになります。

Code 1のときのメモリアクセス
図2 Code 1のときのメモリアクセス
Code 2のときのメモリアクセス
図3 Code 2のときのメモリアクセス

サンプル1と同様にメモリアクセスのジャンプの多いCode 1Code 2より計算時間が長いです。サイズnが大きくなるにつれてCode 1Code 2との計算時間の比率が増えていき、Code 1Code 2より80倍以上遅くなることもあります。



はじめに

演算数を減らす

メモリジャンプを減らす

高性能のアルゴリズム

その他



3 4 9 0 3 0