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

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

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

言語の変更:   FORTRAN版   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
 
/*
  Program to measure time to carry out matrix-vector product.
  c(i) = sum(a(j,i)*b(j),j=1,...,n) for i=1,...,n
 */

#include <stdio.h>
#include <time.h>
int main()
{
  int i, j;
  printf("Matrix size    Elapsed time [sec] \n");

  for (int k=1; k<=30; k++) {
    // Matrix size
    int n = k*1000;

    // Allocation
    int **a, *b, *c;
    a = new int*[n];
    for (i=0; i<n; i++)
      a[i] = new int[n];
    b = new int[n];
    c = new int[n];

    // Initialization
    for (j=0; j<n; j++)
      for (i=0; i<n; i++)
	a[j][i] = 1;
    for (i=0; i<n; i++)
      b[i] = 1;

    // Start time
    clock_t time0 = clock();

    // Main calculation: matrix-vector product
    for (i=0; i<n; i++) {
      c[i] = 0;
      for (j=0; j<n; j++) {
	c[i] += a[j][i] * b[j];
      }
    }

    // Finish time
    clock_t time1 = clock();

    // Output time
    double eTime = (double)(time1-time0)/CLOCKS_PER_SEC;
    printf("%11d %15.7f\n",n, eTime);

    // Deallocation
    for (int i=0; i<n; i++)
      delete[] a[i];
    delete[] a;
    delete[] b;
    delete[] c;
  }
  return 0;
}

    
 
/*
  Program to measure time to carry out matrix-vector product.
  c(i) = sum(a(i,j)*b(j),j=1,...,n) for i=1,...,n
 */

#include <stdio.h>
#include <time.h>
int main()
{
  int i, j;
  printf("Matrix size    Elapsed time [sec] \n");

  for (int k=1; k<=30; k++) {
    // Matrix size
    int n = k*1000;

    // Allocation
    int **a, *b, *c;
    a = new int*[n];
    for (i=0; i<n; i++)
      a[i] = new int[n];
    b = new int[n];
    c = new int[n];

    // Initialization
    for (j=0; j<n; j++)
      for (i=0; i<n; i++)
	a[j][i] = 1;
    for (i=0; i<n; i++)
      b[i] = 1;

    // Start time
    clock_t time0 = clock();

    // Main calculation: matrix-vector product
    for (i=0; i<n; i++) {
      c[i] = 0;
      for (j=0; j<n; j++) {
	c[i] += a[i][j] * b[j];
      }
    }

    // Finish time
    clock_t time1 = clock();

    // Output time
    double eTime = (double)(time1-time0)/CLOCKS_PER_SEC;
    printf("%11d %15.7f\n",n, eTime);

    // Deallocation
    for (int i=0; i<n; i++)
      delete[] a[i];
    delete[] a;
    delete[] b;
    delete[] c;
  }
  return 0;
}

    

■ 計算時間の測定結果

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

測定時間
図1 測定時間

■ 考察

C言語では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より50倍以上遅くなることもあります。



はじめに

演算数を減らす

メモリジャンプを減らす

高性能のアルゴリズム

その他



4 8 2 6 0 7