🚀

AVX2で書いた配列の足し算

2024/04/21に公開

GCC

avx2

tech

はじめに

昔書いた簡単なAVX2を使用したプログラムを見直して、何をしているかを思い出しながら実際に動かしてみたいと思います。

AVX2とは？

AVX2とは「Intel Advanced Vector Extensions 2」の略で、Intel系のCPUで使われる拡張命令セットです。Vactorとついているとおりベクトル演算を可能にしています。
AVX等の拡張命令のデータシートは以下にあります。

余談:自分のマシンがAVX2に対応しているか確認する方法

まずは、自分のマシンのCPUがAVX2に対応しているか確認する方法をいくつか説明します。

Linux の場合は、lscpuなどでCPUの情報を取得すればわかります。

$ lscpu |grep -i avx2
Flags:                                fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf tsc_known_freq pni pclmulqdq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefetch cpuid_fault epb ssbd ibrs ibpb stibp ibrs_enhanced tpr_shadow vnmi flexpriority ept vpid ept_ad fsgsbase tsc_adjust bmi1 avx2 smep bmi2 erms invpcid rdseed adx smap clflushopt clwb intel_pt sha_ni xsaveopt xsavec xgetbv1 xsaves split_lock_detect avx_vnni dtherm ida arat pln pts hwp hwp_notify hwp_act_window hwp_epp hwp_pkg_req hfi umip pku ospke waitpkg gfni vaes vpclmulqdq tme rdpid movdiri movdir64b fsrm md_clear serialize pconfig arch_lbr ibt flush_l1d arch_capabilities

この方法以外にも以下のページのようにすれば良いでしょう。

私のマシンはCPUが13th Gen Intel(R) Core(TM) i5-13600Kなのでありました。

今回の演算

今回のプログラムは単純なベクトルの足し算をします。実行時間の計測にはOpenMPのomp_get_wtime()関数を使用して簡単に取得できるようにします。

プログラムの全体像

プログラムは以下の通りです。

avx2test.c

#include<stdio.h>
#include<immintrin.h>
#include<omp.h>

#define N 5000000
int aa[N*8],bb[N*8],cc[N*8];
__m256i a[N],b[N],c[N];
int main(int argc,char *argv[])
{
        int i,n;
        int *A,*B,*C;
        for(n=0;n<N;n++)
        {
                for(i=0;i<8;i++)
                {
                        aa[i+n*8]=rand()%128;
                        bb[i+n*8]=rand()%128;
                }
                a[n]=_mm256_set_epi32(aa[7+n*8],aa[6+n*8],aa[5+n*8],aa[4+n*8],aa[3+n*8],aa[2+n*8],aa[1+n*8],aa[0+n*8]);
                b[n]=_mm256_set_epi32(bb[7+n*8],bb[6+n*8],bb[5+n*8],bb[4+n*8],bb[3+n*8],bb[2+n*8],bb[1+n*8],bb[0+n*8]);
        }
        double A_wtime=omp_get_wtime();
        for(n=0;n<N;n++)
                for(i=0;i<8;i++)
                        cc[i+n*8]=aa[i+n*8]+bb[i+n*8];
        A_wtime=omp_get_wtime()-A_wtime;
        double V_wtime=omp_get_wtime();
        for(n=0;n<N;n++)
                c[n]=_mm256_add_epi32(a[n],b[n]);
        V_wtime=omp_get_wtime()-V_wtime;
        for(n=0;n<N;n++)
        {
                A=(int*)&a[n];
                B=(int*)&b[n];
                C=(int*)&c[n];
                for(i=0;i<8;i++)
                        if(C[i]!=cc[i+n*8])
                                printf("(false)n=%5d,%3d,%3d=%3d+%3d\n",i+(8*n),cc[i+n*8],C[i],A[i],B[i]);
        }
        printf("Array %lf,Vector %lf\n",A_wtime,V_wtime);
        return 0;
}

では、一個ずつ解説していきましょう。

ヘッダファイルの読み込み

ヘッダファイルはstdio.hの他にOpenMPの関数を使うのでomp.h、AVX2の命令を使うためにはimmintrin.hが必要になります。

#include<stdio.h>
#include<immintrin.h>
#include<omp.h>

繰りかえしの設定

今回は5百万回#define N 5000000実行して差が出るか確認します。

配列の用意

今回はC=A+Bと計算させるのでベクトル演算用にa,b,cと配列演算用にaa,bb,ccを用意しました。また今回は256ビットを同時に使うので以下のようにデータを確保します。

int aa[N*8],bb[N*8],cc[N*8];
__m256i a[N],b[N],c[N];

このときのintのサイズは4バイトだったので、256ビットのベクトルだと8個のデータが一度に計算できるはずです。なのでaa,bb,ccはN*8のサイズを取ります。

データの生成

データはランダムに生成したものの mod 128を取ったものを使います。

for(n=0;n<N;n++)
{
  for(i=0;i<8;i++)
  {
    aa[i+n*8]=rand()%128;
    bb[i+n*8]=rand()%128;
  }
}

こんな変な配列のインデックスを指示しいてるのは特に意味はありません、8個×5,000,000回まわしてデータを格納しているに過ぎません。

データの格納

ではこのa,bにデータをAVX2用に格納していきましょう。
データシートによると、
__m256i _mm256_set_epi32 (int e7, int e6, int e5, int e4, int e3, int e2, int e1, int e0)と言っているので、データがe0,e1,...,e6,e7と並んでいる配列aaならば、
a=_mm256_set_epi32(e7,e6,...,e1,e0);と言うように格納しなさいと言っています。
それがこのようになります。bも同じです。

for(n=0;n<N;n++)
{
  for(i=0;i<8;i++)
  {
    aa[i+n*8]=rand()%128;
    bb[i+n*8]=rand()%128;
  }
　　a[n]=_mm256_set_epi32(aa[7+n*8],aa[6+n*8],aa[5+n*8],aa[4+n*8],aa[3+n*8],aa[2+n*8],aa[1+n*8],aa[0+n*8]);
　　b[n]=_mm256_set_epi32(bb[7+n*8],bb[6+n*8],bb[5+n*8],bb[4+n*8],bb[3+n*8],bb[2+n*8],bb[1+n*8],bb[0+n*8]);
}

計算①:配列での計算

配列の計算は通常通りやります。最初とおわりにomp_get_wtime()関数を入れて実行時間を計測します。

double A_wtime=omp_get_wtime();
for(n=0;n<N;n++)
  for(i=0;i<8;i++)
    cc[i+n*8]=aa[i+n*8]+bb[i+n*8];
A_wtime=omp_get_wtime()-A_wtime;

計算②:ベクトル演算

続いてベクトル演算をします。これも最初とおわりにomp_get_wtime()関数を入れて実行時間を計測します。

double V_wtime=omp_get_wtime();
for(n=0;n<N;n++)
  c[n]=_mm256_add_epi32(a[n],b[n]);
V_wtime=omp_get_wtime()-V_wtime;

検算

結果が合っているか比較して間違っている場合にはどこが間違っているか表示するようにします。

for(n=0;n<N;n++)
{
  A=(int*)&a[n];
  B=(int*)&b[n];
  C=(int*)&c[n];
  for(i=0;i<8;i++)
    if(C[i]!=cc[i+n*8])
      printf("(false)n=%5d,%3d,%3d=%3d+%3d\n",i+(8*n),cc[i+n*8],C[i],A[i],B[i]);
}

結果表示

実行時間がどれだけかかったか最後に表示します。

printf("Array %lf,Vector %lf\n",A_wtime,V_wtime);

コンパイル&実行

では実際にコンパイルして比較してみましょう。

$ gcc -fopenmp -mavx2 -o avx2test avx2test.c
$ ./avx2test
Array 0.069858,Vector 0.025931

効果は出ているみたいですね。
実行順序でキャッシュが効いているといけないので逆でもやりましたが、変わりませんでした。

余談:実行時間が拮抗する場合がある

しかし、こうすると実行時間に差はありませんでした。

$ gcc -O2 -fopenmp -mavx2 -o avx2test avx2test.c
$ ./avx2test
Array 0.029258,Vector 0.029575

つまりGCCの最適化は簡単なベクトル演算程度なら凌駕するということでしょうか？最適化についての詳細は今後の課題です。

おわりに

今回はベクトル演算の初歩である配列をベクトルに変換して足し算するということを書いてみました。今後様々な計算をしていきたいと思います。