BLASを用いた行列積の高速化

Question

前提・実現したいこと

Netlibのサイト上のBLASのdgemm関数ををC言語から呼び出して行列積を計算し、性能比較を行ったのですが、本来BLASを用いたものの方が性能が良くなるはずですが、自分の書いたコードの方が良い性能を出してしまいました。「BLASを用いたコード」の方の誤りを指摘していただきたいです。

発生している問題・エラーメッセージ

以下の表の通り、自分で最適化したコードの方が性能が良くなってしまいます。
性能計算には、行列サイズN、計測時間t[s]として
MFlops = 2n^3 / t
で定義されるMFlopsという値を用いました。

行列サイズ	自分で書いたコードの性能	BLASの性能
1000	4266.67	2000
1100	4259.2	1936
1200	3949.71	2067.140187
1300	4260.85	2008.685714
1400	4231.71	1973.213483
1500	4500	1963.636364
1600	4096	2032.124031
1700	4249.08	1894.168675
1800	4496.96	1852.347395
1900	4456.61	1817.706004
2000	4551.11	1865.209472
2100	4374.2	1863.849057
2200	4340.59	1831.913978
2300	4338.09	1830.054054
2400	4243.34	1938.085433
2500	4366.81	1913.875598
2600	4285.2	1973.445614
2700	4321.48	2036.721099
2800	4383.55	1975.988748
2900	4421.8	2007.58328
3000	4330.83	2003.478261

該当のソースコード

自分で最適化したコード

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 3000 // Max matrix size
#define M 50 // Block size

double A[N][N], B[N][N], C[N][N];

void mm_ikj_opt(int, int, int);

int main()
{
  int i, j, n;
  double t, MFlops;
  clock_t tic, toc;

  int index;
  printf("N,t,MFlops\n");

  for(index = 1000; index <= 3000; index += 100) { // 行列サイズ1000~3000まで100おきにデータを取る
    n = index;
    srand(2021);
    for(i = 0; i < n; i++) {
      for(j = 0; j < n; j++) {
        A[i][j] = rand();
        B[i][j] = rand();
        C[i][j] = 0.0;
      }
    }

    tic = clock();
    mm_ikj_opt(n, n, n);
    toc = clock();
    t = (double)(toc - tic) / CLOCKS_PER_SEC;
    MFlops = (2.0 * n) * n * n / t / 1e+6; // 性能
    printf("%d, %.2f, %.2f\n", n, t, MFlops);
  }
  return 0;
}

void mm_ikj_opt(int m, int p, int n)
{
  int ii, jj, kk;
  int i, k, j;

  // ループアンローリング、ブロック化などで最適化
  for(ii = 0; ii < m; ii += M) {
    for(kk = 0; kk < p; kk += M) {
      for(jj = 0; jj < n; jj += M) {
        for(i = ii; i < ii + M; i++) {
          for(k = kk; k < kk + M; k++) {
            for(j = jj; j < jj + M - 4; j += 4) { // loop unrolling
              C[i][j] += A[i][k] * B[k][j];
              C[i][j + 1] += A[i][k] * B[k][j + 1];
              C[i][j + 2] += A[i][k] * B[k][j + 2];
              C[i][j + 3] += A[i][k] * B[k][j + 3];
            }
          }
        }
      }
    }
  }
}

BLASを用いたコード

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 3000 // matrix size
double A[N * N], B[N * N], C[N * N];

int main()
{
  int i, j, n;
  double t, MFlops;
  clock_t tic, toc;

  int index;
  printf("N,t,MFlops\n");

  for(index = 1000; index <= 3000; index += 100) {
    n = index;
    srand(2021);
    double alpha = 1.0, beta = 0.0;
    for(i = 0; i < n; i++) {
      for(j = 0; j < n; j++) {
        int id1 = i * n + j;
        int id2 = i + j * n;
        A[id1] = rand();
        B[id2] = rand();
        C[id1] = 0.0;
      }
    }

    tic = clock();
    dgemm_("T", "N", &n, &n, &n, &alpha, A, &n, B, &n, &beta, C, &n); // dgemmで行列積を計算
    toc = clock();
    t = (double)(toc - tic) / CLOCKS_PER_SEC;
    MFlops = (2.0 * n) * n * n / t / 1e+6;
    printf("%d, %f, %f\n", n, t, MFlops);

  }
  return 0;
}

試したこと

FORTRANが列方向アクセスなのに対してCは行方向アクセスなので、行方向に計算を行っていく行列Aの方を値を代入する段階で転置した状態で代入し、dgemm関数内でも行列Aを転置させて計算することで速くなると思ったが、ならなかった。

補足情報（FW/ツールのバージョンなど）

OS: Windows 11
CPU: AMD Ryzen 7 3700U with Radeon Vega Mobile Gfx 2.30 GHz
コンパイラ: gcc (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0

コンパイル・実行方法

$ gcc -O2 hoge.c -lblas -o hoge
$ ./hoge > hoge.csv

Comments

[jbpb0]

> 本来BLASを用いたものの方が性能が良くなるはず

https://note.com/ymzo76/n/ncca5e176c15e#0TzlF
の「性能比較」が参考になるかもしれません

[episteme]

> 本来BLASを用いたものの方が性能が良くなるはずです

その根拠は?

[PC_breakman]

FORTRANで書かれたBLASにCからアクセスしようとして遅くなっていたのかもしれません。CBLASを使ったら最適化の2倍以上の性能を出せました。ありがとうございます。

Answer 1

CBLASを用いたら2倍以上の性能が出せました。