6要素AVX

Question

で倍精度ベクトル行列ベクトル乗算Iは倍精度で以下の操作を実行する必要があります

数字は値がメモリに格納される方法を表します。私はこれをAVXで実装したいと思います。 [QK]の欄を最初に8つの要素にパディングしてから、[x]と[QK]でマトリックスベクトル乗算を実行した後、ドットプロダクトを実行するのが最善でしょうか？

EDIT：[OK]を、私は以下に示すようにパディングベクターをFLOAT32-ビットバージョンを実装することを決定した：

// Perform matrix vector multiplication of QK*x    
// Load first four columns QK into 4 ymm registers  
ymm0 = _mm256_load_ps((float *)(QK));  
ymm1 = _mm256_load_ps((float *)(QK+8));  
ymm2 = _mm256_load_ps((float *)(QK+16));  
ymm3 = _mm256_load_ps((float *)(QK+24)); 

// Load first four values of x 
ymm4 = _mm256_broadcast_ss((float *)(x)); 
ymm5 = _mm256_broadcast_ss((float *)(x+1)); 
ymm6 = _mm256_broadcast_ss((float *)(x+2)); 
ymm7 = _mm256_broadcast_ss((float *)(x+3)); 

// Multiply in place - frees up ymm4,ymm5,ymm6,ymm7 
ymm0 = _mm256_mul_ps(ymm0, ymm4); 
ymm1 = _mm256_mul_ps(ymm1, ymm5); 
ymm2 = _mm256_mul_ps(ymm2, ymm6); 
ymm3 = _mm256_mul_ps(ymm3, ymm7); 

// Add together, this frees up ymm1,ymm2,and ymm3 
ymm0 = _mm256_add_ps(ymm0, ymm1); 
ymm2 = _mm256_add_ps(ymm2, ymm3); 
ymm0 = _mm256_add_ps(ymm0, ymm2); 

// Load next last columns of QK 
ymm1 = _mm256_load_ps((float *)(QK+32));  
ymm2 = _mm256_load_ps((float *)(QK+40)); 

// Load last two values of x 
ymm6 = _mm256_broadcast_ss((float *)(x+4)); 
ymm7 = _mm256_broadcast_ss((float *)(x+5)); 

// Multiply in place 
ymm1 = _mm256_mul_ps(ymm1, ymm6); 
ymm2 = _mm256_mul_ps(ymm2, ymm7); 

// Add together, this frees up every register except for ymm0 
ymm0 = _mm256_add_ps(ymm0, ymm1); 
ymm0 = _mm256_add_ps(ymm0, ymm2); 

// Answer stored in ymm0 and ymm1 
// Calculate dot product of y*(QK*x) 
// Load x 
ymm1 = _mm256_load_ps((float *)(y)); 

// Do dotproduct by using horizontal multiply followed by horizontal add 
// Multiply in place 
ymm0 = _mm256_mul_ps(ymm0, ymm1); 

// Do horizontal sum 
__m128 xmm1 = _mm256_extractf128_ps(ymm0, 1); 
__m128 xmm2 = _mm256_extractf128_ps(ymm0, 0); 
xmm2 = _mm_add_ps(xmm1, xmm2); 
xmm1 = _mm_movehl_ps(xmm1, xmm2); 
xmm2 = _mm_add_ps(xmm1, xmm2); 
xmm1 = _mm_shuffle_ps(xmm2, xmm2, 1); 
xmm2 = _mm_add_ss(xmm1, xmm2); 
ans[0] = _mm_cvtss_f32(xmm2); 

現状では、それは約3倍の速さよりも実行されます。

ans[0] = (QK[0]*x[0]+QK[8]*x[1]+QK[16]*x[2]+QK[24]*x[3]+QK[32]*x[4]+QK[40]*x[5])*y[0]+ 
     (QK[1]*x[0]+QK[9]*x[1]+QK[17]*x[2]+QK[25]*x[3]+QK[33]*x[4]+QK[41]*x[5])*y[1]+ 
     (QK[2]*x[0]+QK[10]*x[1]+QK[18]*x[2]+QK[26]*x[3]+QK[34]*x[4]+QK[42]*x[5])*y[2]+ 
     (QK[3]*x[0]+QK[11]*x[1]+QK[19]*x[2]+QK[27]*x[3]+QK[35]*x[4]+QK[43]*x[5])*y[3]+ 
     (QK[4]*x[0]+QK[12]*x[1]+QK[20]*x[2]+QK[28]*x[3]+QK[36]*x[4]+QK[44]*x[5])*y[4]+ 
     (QK[5]*x[0]+QK[13]*x[1]+QK[21]*x[2]+QK[29]*x[3]+QK[37]*x[4]+QK[45]*x[5])*y[5]; 

繰り返し回数が500回の場合、標準Cバージョンは約9秒で実行され、単精度AVXバージョンは約3.5秒で実行されます。私が最後に水平の合計をコメントアウトすると、それは約0.5秒で実行されます。水平合計は実際にパフォーマンスを殺します...

Answer 1

これを効率的に行うためのコードを作成しました。私は、AVXをシングルスレッド用に使用すると、ほぼ4倍のスピードを上げます（最高の倍率でAVXに期待できます）。ここには、6x6マトリックス上の2000、32000、および4000000のxとyの6成分ベクトルで実行される結果があります。これらはおおよそ私のシステムのL2、L3、および>> L3キャッシュサイズに対応しています（各ベクトルは48バイトを使用します）。

編集：フロートでこれを行うには、最後にテキスト（およびコード）を追加しました。単一のスレッドでAVXを使用するとスピードアップが8倍近くになります。

i5-3317U (2 core ivy bridge) 
compiled with: g++ m6.cpp -o m6 -O3 -mavx -fopenmp 
nvec 2000, repeat 10000, 1 thread : time scalar/SIMD 3.95 
nvec 32000, repeat 1000, 1 thread : time scalar/SIMD 3.53 
nvec 4000000, repeat 10, 1 thread : time scalar/SIMD 3.28 
1 thread for both the SIMD and non-SIMD functions 

nvec 2000, repeat 10000, 2 threads: time scalar/SIMD 5.96 
nvec 32000, repeat 1000, 2 threads: time scalar/SIMD 5.88 
nvec 4000000, repeat 10, 2 threads: time scalar/SIMD 4.52 
2 threads on the SIMD function vs. 1 thread on the non-SIMD function 

compiled with g++ m6.cpp -o m6 -O3 -msse4.2 -fopenmp 
nvec 2000, repeat 10000, 1 thread: time scalar/SIMD 2.15 
nvec 32000, repeat 1000, 1 thread: time scalar/SIMD 2.06 
nvec 4000000, repeat 10, 1 thread: time scalar/SIMD 2.00 

基本アルゴリズムは、6x6行列ではなく、xおよびyベクトル上でSIMDを実行します。 1つの重要な点は、xとyのベクトルは配列の構造体のブロックの配列でなければならないということです。これは、配列構造体の配列（AoSoA）または配列のハイブリッド構造体とも呼ばれます。ここで詳しく読むことができます。「携帯型SIMDプログラミングのためのC言語の拡張」 http://www.sci.utah.edu/~wald/Publications/2012///ppopp/download//vecimp.pdf

以下はコードです。 aos2aosoa関数は、6つのコンポーネントベクトルの配列をSoAの配列に変換します。 SIMDの利点を最大限に活用したい場合、アプリケーションはこの形式のベクトルを生成する必要があります（変換は行いません）。このコードでは、Agner Fogのベクトルクラスhttp://www.agner.org/optimize/#vectorclassを使用しています。それはちょうどいくつかのヘッダーファイルです。このコードはSSEでも動作します（ただし、ブーストは予想どおり約2倍です）。

xとyのベクトルの配列と結果は4の倍数でなければなりません。ただし、ベクトルの数はありません。コンパイラのcommmandラインオプションでAVX

コード：

のVisual Studio PUT /アーチでこの

g++ m6.cpp -o m6 -O3 -mavx -fopenmp //system with AVX 
g++ m6.cpp -o m6 -O3 -msse4.2 -fopenmp //system with SSE 

のようなコンパイルここで

#include <stdio.h> 
#include <omp.h> 
#include "vectorclass.h" 
#include <stdlib.h> 

double prod_scalar(double *x, double *M, double *y) { 
    double sum = 0.0f; 
    for(int i=0; i<6; i++) { 
     for(int j=0; j<6; j++) { 
      sum += x[i]*M[i*6 + j]*y[j]; 
     } 
    } 
    return sum; 
} 

void prod_block4(double *x, double *M, double *y, double *result) { 
    Vec4d sum4 = 0.0f; 
    for(int i=0; i<6; i++) { 
     Vec4d x4 = Vec4d().load(&x[4*i]); 
     for(int j=0; j<6; j++) { 
      Vec4d y4 = Vec4d().load(&y[4*j]); 
      sum4 += x4*M[i*6 + j]*y4; 
     } 
    } 
    sum4.store(result); 
} 

void prod_block4_unroll2(double *x, double *M, double *y, double *result) { 
    Vec4d sum4_1 = 0.0f; 
    Vec4d sum4_2 = 0.0f; 
    Vec4d yrow[6]; 
    for(int i=0; i<6; i++) { 
     yrow[i] = Vec4d().load(&y[4*i]); 
    } 
    for(int i=0; i<6; i++) { 
     Vec4d x4 = Vec4d().load(&x[4*i]); 
     for(int j=0; j<6; j+=2) { 
      sum4_1 += x4*M[i*6 + j]*yrow[j]; 
      sum4_2 += x4*M[i*6 + j+1]*yrow[j+1]; 
     } 
    } 
    sum4_1 += sum4_2; 
    sum4_1.store(result); 
} 
void loop_scalar(double *x, double *M, double *y, double *result, const int nvec) { 
// #pragma omp parallel for 
    for(int i=0; i<nvec; i++) { 
     result[i] = prod_scalar(&x[6*i], M, &y[6*i]); 
    } 
} 

void loop_block4(double *x, double *M, double *y, double *result, const int nvec) { 
// #pragma omp parallel for 
    for(int i=0; i<(nvec/4); i++) { 
//  prod_block4(&x[i*24], M, &y[i*24], &result[4*i]); 
     prod_block4_unroll2(&x[i*24], M, &y[i*24], &result[4*i]); 
    } 
} 


void aos2soa(double *in, double *out) { 
    int cnt = 0; 
    for(int i=0; i<6; i++) { 
     for(int j=0; j<4; j++) { 
      out[cnt++] = in[6*j + i]; 
     } 
    } 
} 

void aos2aosoa(double *in, double *out, const int nvec) { 
    for(int i=0; i<(nvec/4); i++) { 
     aos2soa(&in[i*24], &out[i*24]); 
    } 
} 

double compare_results(double *r1, double * r2, const int nvec) { 
    double out = 0.0f; 
    for(int i=0; i<nvec; i++) { 
     out += r1[i] -r2[i]; 
     //if(out!=0) printf("%f %f\n",r1[i], r2[i]); 
    } 
    return out; 
} 

void loop(const int nvec, const int repeat) { 
    double *M = (double*)_mm_malloc(6*6*sizeof(double),32); 
    double *x_aos = (double*)_mm_malloc(nvec*6*sizeof(double),32); 
    double *y_aos = (double*)_mm_malloc(nvec*6*sizeof(double),32); 
    double *x_aosoa = (double*)_mm_malloc(nvec*6*sizeof(double),32); 
    double *y_aosoa = (double*)_mm_malloc(nvec*6*sizeof(double),32); 
    double *results_scalar = (double*)_mm_malloc(nvec*sizeof(double),32); 
    double *results_vector = (double*)_mm_malloc(nvec*sizeof(double),32); 

    for(int i=0; i<6; i++) { 
     for(int j=0; j<6; j++) { 
      M[j*6 + i] = i*6 + j; 
     } 
    } 
    for(int i=0; i<(6*nvec); i++) { 
     double r1 = (double)rand()/RAND_MAX; 
     double r2 = (double)rand()/RAND_MAX; 
     //x_aos[i] = i; 
     x_aos[i] = r1; 
     //y_aos[i] = i; 
     y_aos[i] = r2; 

    } 

    aos2aosoa(x_aos, x_aosoa, nvec); 
    aos2aosoa(y_aos, y_aosoa, nvec); 

    double dtime; 
    dtime = omp_get_wtime(); 
    for(int i=0; i<repeat; i++) { 
     loop_scalar(x_aos, M, y_aos, results_scalar, nvec); 
    } 
    dtime = omp_get_wtime() - dtime; 
    printf("time scalar %f\n", dtime); 

    double dtime_old = dtime; 
    dtime = omp_get_wtime(); 
    for(int i=0; i<repeat; i++) { 
     loop_block4(x_aosoa, M, y_aosoa, results_vector, nvec); 
    } 
    dtime = omp_get_wtime() - dtime; 
    printf("time vector %f\n", dtime); 
    printf("time scalar/vector %f\n", dtime_old/dtime); 

    printf("difference %f\n", compare_results(results_scalar, results_vector, nvec)); 

    _mm_free(M); 
    _mm_free(x_aos); 
    _mm_free(y_aos); 
    _mm_free(x_aosoa); 
    _mm_free(y_aosoa); 
    _mm_free(results_scalar); 
    _mm_free(results_vector); 

} 

int main() { 
    int nveca[3]; 
    nveca[0] = 2000; // 2000*2*6*8 = 192kb //L2 
    nveca[1] = 32000; // 32000*2*6*8 = 3Mb //L3 
    nveca[2] = 4*1000000; //366Mb 
    int nrepeat[3]; 
    nrepeat[0] = 10000; 
    nrepeat[1] = 1000; 
    nrepeat[2] = 10; 
    for(int i=0; i<3; i++) { 
     printf("nvec %d, repeat %d\n", nveca[i], nrepeat[i]); 
     loop(nveca[i], nrepeat[i]); 
     printf("\n"); 
    } 
} 

フロートための結果があります。ブーストはおよそここ

nvec 2000, repeat 10000, 1 thread: time scalar/SIMD 7.86 
nvec 32000, repeat 1000, 1 thread: time scalar/SIMD 7.63 
nvec 5000000, repeat 100, 1 thread: time scalar/SIMD 6.63 

8倍速フロートのコードは、あなたがベクトルの行列要素1のためのYMMレジスタの8 ---> 6を使用することができますWindows 7の代わりに、二重

#include <stdio.h> 
#include <omp.h> 
#include "vectorclass.h" 
#include <stdlib.h> 

#define ROUND_DOWN(x, s) ((x) & ~((s)-1)) 

float prod_scalar(float *x, float *M, float *y) { 
    float sum = 0.0f; 
    for(int i=0; i<6; i++) { 
     for(int j=0; j<6; j++) { 
      sum += x[i]*M[i*6 + j]*y[j]; 
     } 
    } 
    return sum; 
} 

float prod_scalar_unroll2(float *x, float *M, float *y) { 
    float sum_1 = 0.0f; 
    float sum_2 = 0.0f; 
    for(int i=0; i<6; i++) { 
     for(int j=0; j<6; j+=2) { 
      sum_1 += x[i]*M[i*6 + j]*y[j]; 
      sum_2 += x[i]*M[i*6 + j+1]*y[j+1]; 
     } 
    } 
    return sum_1 + sum_2; 
} 

void prod_block4(float *x, float *M, float *y, float *result) { 
    Vec8f sum4 = 0.0f; 
    for(int i=0; i<6; i++) { 
     Vec8f x4 = Vec8f().load(&x[8*i]); 
     for(int j=0; j<6; j++) { 
      Vec8f y4 = Vec8f().load(&y[8*j]); 
      sum4 += x4*M[i*6 + j]*y4; 
     } 
    } 
    sum4.store(result); 
} 

void prod_block4_unroll2(float *x, float *M, float *y, float *result) { 
    Vec8f sum4_1 = 0.0f; 
    Vec8f sum4_2 = 0.0f; 
    Vec8f yrow[6]; 
    for(int i=0; i<6; i++) { 
     yrow[i] = Vec8f().load(&y[8*i]); 
    } 
    for(int i=0; i<6; i++) { 
     Vec8f x4 = Vec8f().load(&x[8*i]); 
     for(int j=0; j<6; j+=2) { 
      sum4_1 += x4*M[i*6 + j]*yrow[j]; 
      sum4_2 += x4*M[i*6 + j+1]*yrow[j+1]; 
     } 
    } 
    sum4_1 += sum4_2; 
    sum4_1.store(result); 
} 

void loop_scalar(float *x, float *M, float *y, float *result, const int nvec) { 
// #pragma omp parallel for 
    for(int i=0; i<nvec; i++) { 
     result[i] = prod_scalar(&x[6*i], M, &y[6*i]); 
     //result[i] = prod_scalar_unroll2(&x[6*i], M, &y[6*i]); 
    } 
} 

void loop_SIMD(float *x, float *M, float *y, float *result, const int nvec) { 
    //#pragma omp parallel for schedule(static,256) 
    //#pragma omp parallel for schedule(static) 
    const int N = nvec/8; 
    //printf("chuck %d\n", N/2); 
    //omp_set_num_threads(2); 

    //#pragma omp parallel 
    { 
     //int nthreads = omp_get_num_threads(); 
     //int ithread = omp_get_thread_num(); 

     //int start = (ithread*N)/nthreads; 
     //int end = ((ithread+1)*N)/nthreads; 
     //printf("ithread, start %d, end %d, chunk %d\n",start, end, end-start); 
     //#pragma omp for 
     for(int i=0; i<(nvec/8); i++) { 
     //for(int i=start; i<end; i++) { 
    //  prod_block4(&x[i*24], M, &y[i*24], &result[4*i]); 
     prod_block4_unroll2(&x[i*48], M, &y[i*48], &result[8*i]); 
     } 
    } 
} 

void aos2soa(float *in, float *out) { 
    int cnt = 0; 
    for(int i=0; i<6; i++) { 
     for(int j=0; j<8; j++) { 
      out[cnt++] = in[6*j + i]; 
     } 
    } 
} 

void aos2aosoa(float *in, float *out, const int nvec) { 
    for(int i=0; i<(nvec/8); i++) { 
     aos2soa(&in[i*48], &out[i*48]); 
    } 
} 

float compare_results(float *r1, float * r2, const int nvec) { 
    float out = 0.0f; 
    for(int i=0; i<nvec; i++) { 
     out += r1[i] -r2[i]; 
     //if(out!=0) printf("%f %f\n",r1[i], r2[i]); 
    } 
    return out; 
} 

void loop(const int nvec, const int repeat) { 
    float *M = (float*)_mm_malloc(6*6*sizeof(float),64); 
    float *x_aos = (float*)_mm_malloc(nvec*6*sizeof(float),64); 
    float *y_aos = (float*)_mm_malloc(nvec*6*sizeof(float),64); 
    float *x_aosoa = (float*)_mm_malloc(nvec*6*sizeof(float),64); 
    float *y_aosoa = (float*)_mm_malloc(nvec*6*sizeof(float),64); 
    float *results_scalar = (float*)_mm_malloc(nvec*sizeof(float),64); 
    float *results_SIMD = (float*)_mm_malloc(nvec*sizeof(float),64); 

    for(int i=0; i<6; i++) { 
     for(int j=0; j<6; j++) { 
      M[j*6 + i] = i*6 + j; 
     } 
    } 
    for(int i=0; i<(6*nvec); i++) { 
     float r1 = (float)rand()/RAND_MAX; 
     float r2 = (float)rand()/RAND_MAX; 
     //x_aos[i] = i; 
     x_aos[i] = r1; 
     //y_aos[i] = i; 
     y_aos[i] = r2; 

    } 

    aos2aosoa(x_aos, x_aosoa, nvec); 
    aos2aosoa(y_aos, y_aosoa, nvec); 

    float dtime; 
    dtime = omp_get_wtime(); 
    for(int i=0; i<repeat; i++) { 
     loop_scalar(x_aos, M, y_aos, results_scalar, nvec); 
    } 
    dtime = omp_get_wtime() - dtime; 
    printf("time scalar %f\n", dtime); 

    float dtime_old = dtime; 
    dtime = omp_get_wtime(); 
    for(int i=0; i<repeat; i++) { 
     loop_SIMD(x_aosoa, M, y_aosoa, results_SIMD, nvec); 
    } 
    dtime = omp_get_wtime() - dtime; 
    printf("time SIMD %f\n", dtime); 
    printf("time scalar/SIMD %f\n", dtime_old/dtime); 

    printf("difference %f\n", compare_results(results_scalar, results_SIMD, nvec)); 

    _mm_free(M); 
    _mm_free(x_aos); 
    _mm_free(y_aos); 
    _mm_free(x_aosoa); 
    _mm_free(y_aosoa); 
    _mm_free(results_scalar); 
    _mm_free(results_SIMD); 

} 

int main() { 
    int nveca[3]; 
    nveca[0] = 2000; // 2000*2*6*8 = 192kb //L2 
    nveca[1] = 32000; // 32000*2*6*8 = 3Mb //L3 
    nveca[2] = 5*1000000; //366Mb 
    int nrepeat[3]; 
    nrepeat[0] = 10000; 
    nrepeat[1] = 1000; 
    nrepeat[2] = 100; 
    for(int i=0; i<3; i++) { 
     printf("nvec %d, repeat %d\n", nveca[i], nrepeat[i]); 
     loop(nveca[i], nrepeat[i]); 
     printf("\n"); 
    } 
}