大規模な密行列の乗算にループ タイル/ループ ブロッキングを効果的に使用する方法を誰かが教えてくれるかどうか疑問に思っていました。1000x1000 行列でC = ABを実行しています。ウィキペディアのループ タイリングの例に従いましたが、タイリングを使用しない場合よりも悪い結果が得られます。
http://en.wikipedia.org/wiki/Loop_tiling
以下にいくつかのコードを提供しました。単純な方法は、キャッシュ ミスが原因で非常に遅くなります。transpose メソッドは、バッファ内にBの転置を作成します。この方法は最速の結果をもたらします (行列の乗算は O(n^3) として行われ、転置は O(n^2) として行われるため、転置は少なくとも 1000 倍速くなります)。ブロッキングなしの wiki メソッドも高速であり、バッファは必要ありません。ブロック方法は低速です。ブロッキングのもう 1 つの問題は、ブロックを数回更新する必要があることです。注意しないと競合状態が発生する可能性があるため、これはスレッド化/OpenMP の課題です。
転置法を変更して AVX を使用すると、Eigen よりも速く結果が得られることを指摘しておく必要があります。ただし、SSE での結果は Eigen よりも少し遅いので、キャッシュをより適切に使用できると思います。
編集：私は自分が何をしたいのか考えていると思います。これは、1969 年の Cannon のアルゴリズムに由来します。
http://en.wikipedia.org/wiki/Matrix_multiplication#Communication-avoiding_and_distributed_algorithms
ブロック行列の乗算を行い、各スレッドでAとBではなくCの部分行列を処理する必要があります。たとえば、行列を 4 つのブロックに分割するとします。私はするだろう：
+-+ +-+ +-+ +-+ +-+ +-+
| | | | | |
| C1 C2 | | A1 A2 | | B1 B2 |
| | = | | x | |
| C3 C4 | | A3 A4 | | B3 B4 |
| | | | | |
+-+ +-+ +-+ +-+ +-+ +-+
thread 1: C1 = A1B1 + A2B3
thread 2: C2 = A1B2 + A2B4
thread 3: C3 = A3B1 + A4B3
thread 4: C4 = A3B2 + A4B4
これにより、競合状態が解消されます。私はこれについて考えなければならないでしょう。
void matrix_mult_naive(const float*A , const float* B, float* C, const int N, const int M, const int K) {
#pragma omp parallel for
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
float tmp = 0;
for(int l=0; l<M; l++) {
tmp += A[M*i+l]*B[K*l+j];
}
C[K*i + j] = tmp;
}
}
}
void matrix_mult_transpose(const float*A , const float* B, float* C, const int N, const int M, const int K) {
float *B2 = (float*)aligned_malloc(M*K*sizeof(float), 32);
transpose(B, B2, M, K, 1);
#pragma omp parallel for
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
float tmp = 0;
for(int l=0; l<M; l++) {
tmp += A[M*i+l]*B2[M*j+l];
}
C[K*i + j] = tmp;
}
}
aligned_free(B2);
}
void matrix_mult_wiki(const float*A , const float* B, float* C, const int N, const int M, const int K) {
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
C[K*i + j] = 0;
}
}
#pragma omp parallel for
for(int i=0; i<N; i++) {
for(int l=0; l<M; l++) {
float a = A[M*i+l];
for(int j=0; j<K; j++) {
C[K*i + j] += a*B[K*l+j];
}
}
}
}
void matrix_mult_wiki_block(const float*A , const float* B, float* C, const int N, const int M, const int K) {
const int block_size = 8; //I have tried several different block sizes
for(int i=0; i<N; i++) {
for(int j=0; j<K; j++) {
C[K*i + j] = 0;
}
}
for(int l2=0; l2<M; l2+=block_size) {
for(int j2=0; j2<K; j2+=block_size) {
#pragma omp parallel for
for(int i=0; i<N; i++) {
for(int l=l2; l<min(M, l2+block_size); l++) {
for(int j=j2; j<min(K, j2+block_size); j++) {
C[K*i + j] += A[M*i+l]*B[K*l+j];
}
}
}
}
}
}