

C言語初心者です。  
SIMD命令を使って8×8行列を計算するプログラムを書いたところ、"segmentation fault"
が実行の度に起きたり起きなかったりする奇妙なことになってしまいました。  
どこが間違いなのか、どう改善すべきなのか、教えていただけると幸いです。

またswitch文のところは、ポインタとiを使って計算すべき所だと思うのですが、なぜか変数を使うと "segmentation error"
が起きてしまいます。これに関しても理由が知りたいです。

以下がコードです。

[code]

    #include <math.h>
    #include <stdio.h>
    #include <stdlib.h>
    #include <string.h>
    #include <x86intrin.h>
    
    void print_vec(__m256 m)
    {
       printf("{");
       for(int i=0;i<8; i++) printf("%.2f ",m[i]);
       printf("}\n");
    }
    
    void prod_avx(float* a, float * b, float* c)
    {
      float cd[64];
      __m256 b0n, b1n, b2n, b3n, b4n, b5n, b6n, b7n;
      
      b0n = _mm256_load_ps(b);
      b1n = _mm256_load_ps(b + 8);
      b2n = _mm256_load_ps(b + 16);
      b3n = _mm256_load_ps(b + 24);
      b4n = _mm256_load_ps(b + 32);
      b5n = _mm256_load_ps(b + 40);
      b6n = _mm256_load_ps(b + 48);
      b7n = _mm256_load_ps(b + 56);
    
      for (int i = 0; i < 8; i++)
      {
        __m256 ai0, ai1, ai2, ai3 , ai4, ai5, ai6, ai7;
        ai0 = _mm256_set1_ps(*(a + 8 * i));
        ai1 = _mm256_set1_ps(*(a + 8 * i+1));
        ai2 = _mm256_set1_ps(*(a + 8 * i+2));
        ai3 = _mm256_set1_ps(*(a + 8 * i+3));
        ai4 = _mm256_set1_ps(*(a + 8 * i+4));
        ai5 = _mm256_set1_ps(*(a + 8 * i+5));
        ai6 = _mm256_set1_ps(*(a + 8 * i+6));
        ai7 = _mm256_set1_ps(*(a + 8 * i+7));
        
        ai0 = _mm256_mul_ps(ai0, b0n); 
        ai1 = _mm256_mul_ps(ai1, b1n); 
        ai2 = _mm256_mul_ps(ai2, b2n); 
        ai3 = _mm256_mul_ps(ai3, b3n); 
        ai4 = _mm256_mul_ps(ai4, b4n); 
        ai5 = _mm256_mul_ps(ai5, b5n); 
        ai6 = _mm256_mul_ps(ai6, b6n); 
        ai7 = _mm256_mul_ps(ai7, b7n); 
      
        ai0 = _mm256_add_ps(ai0, ai1);
        ai2 = _mm256_add_ps(ai2, ai3);
        ai4 = _mm256_add_ps(ai4, ai5);
        ai6 = _mm256_add_ps(ai6, ai7);
        
        ai0 = _mm256_add_ps(ai0, ai2);
        ai4 = _mm256_add_ps(ai4, ai6);
        
        ai0 = _mm256_add_ps(ai0, ai4);
    
        //print_vec(ai0);
        //printf("i*8 = %d \n",i*8);
        //理由は本当に不明としか言いようがないが、iをポインタ演算に使わずこうしないと-11エラーを起こす
    
       switch(i)
        {
        case 0: 
            _mm256_store_ps(cd+0 , ai0);
            break;
        case 1: 
            _mm256_store_ps(cd+8 , ai0);
            break;
        case 2:  
            _mm256_store_ps(cd+16 , ai0);
            break;
        case 3:
            _mm256_store_ps(cd+24 , ai0);
            break;
        case 4: 
            _mm256_store_ps(cd+32 , ai0);
            break;
        case 5: 
            _mm256_store_ps(cd+40 , ai0);
            break;
        case 6: 
            _mm256_store_ps(cd+48 , ai0);
            break;
        case 7: 
            _mm256_store_ps(cd+56 , ai0);
            break;
            }
             
     }
      memcpy(c, cd, 64 * sizeof(float));
    }
    void print_gyoretu(float *a)
    { 
     for(int i=0 ; i< 8 ;i++)
        {
        printf("{");
         for(int j=0;j < 8 ; j++)
         {
           printf("%.2f ",a[i*8+j]);  
         }
         printf("}\n");
        }
    }
    
    int main() {
       // float *a=(float *)malloc(sizeof(float) *64);
       // float *b=(float *)malloc(sizeof(float) *64);
        float a[64];
        float b[64];
        float cd[64];
        int j=0;
    
        for(int i=0;i<64;i++)
          {
              if(i%8<4 && i<32)
              {
                a[i]=j;
                j++;
              }
              else
              {
                a[i]=0;
              }
          }
        for(int i=0;i<64;i++)
          {
              if(i%8<4 && i<32)
              {
                b[i]=j;
                j++;
              }
              else
              {
                b[i]=0;
              }
          }
                 
        print_gyoretu(a);
        printf("===========================\n");
        print_gyoretu(b);
    
        printf("\n====================\n");
    
        prod_avx(a,b,cd);
        print_gyoretu(cd);
        
    }
    
[/code]


CC BY-SA 4.0
2021-01-05T12:23:32.960
73097
2021-01-09T01:00:06.340
2021-01-05T12:47:15.567
3060
43429
question
c
SIMDプログラミングによる行列積について


macOS上でXcodeのclangで実行すると事象が再現できないため、他にも問題がある可能性がありますが、`_mm256_load_ps`などの命令は引数のポインタが32バイト境界に整列(アラインメント)していることが必須であるため、たまたま32バイト境界からずれた時にCPU例外("segmentation
fault")が発生している可能性が高いと思います。

うまく整列していない可能性のある配列の宣言を以下のように書き換えてみてください。

`prod_avx`内の`cd`の宣言:

[code]

        float *cd = (float *)aligned_alloc(32, 64 * sizeof(float));
    
[/code]

`main`中の`a`, `b`, `cd`の宣言:

[code]

        float *a = (float *)aligned_alloc(32, 64 * sizeof(float));
        float *b = (float *)aligned_alloc(32, 64 * sizeof(float));
        float *cd = (float *)aligned_alloc(32, 64 * sizeof(float));
    
[/code]

(それぞれ、対応する`free`が必要ですが、簡単のために省略してあります。ご自身で適切に補ってください。)

※gccなら拡張アトリビュートの`__attribute__
((aligned(32)))`が使えるかと思ったのですが、環境によっては32バイトのような大きな境界への整列には対応していない可能性があるため、`aligned_alloc`にしました。


CC BY-SA 4.0
2021-01-05T13:52:36.843
73098
2021-01-05T13:52:36.843
13972
73097
answer


gcc には `-fsanitize=alignment` というオプションがあります。

> -fsanitize=alignment
>
> This option enables checking of alignment of pointers when they are
> dereferenced, or when a reference is bound to insufficiently aligned target,
> or when a method or constructor is invoked on insufficiently aligned object.

このオプションを付けて実行してみると `type '__m256', which requires 32 byte alignment`
と表示されます(SIMD 命令として AVX2 を指定しています)。

[code]

    $ lscpu | grep -E '^(Architecture|Model name)'
    Architecture: x86_64
    Model name:   Intel(R) Core(TM) i5-8500T CPU @ 2.10GHz
    $ lsb_release -ir
    Distributor ID: Ubuntu
    Release:    20.10
    $ gcc --version
    gcc (Ubuntu 10.2.0-13ubuntu1) 10.2.0
    
    # using AVX2
    $ gcc -fsanitize=alignment -mavx2 -Wall -Wextra -g matrix_product.c -o matrix_product && ./matrix_product
                              :
    
    /usr/lib/gcc/x86_64-linux-gnu/10/include/avxintrin.h:874:10:
    runtime error:
    load of misaligned address 0x7ffee9609fd0 for type '__m256', which requires 32 byte alignment
    0x7ffee9609fd0: note: pointer points here
     00 00 00 00  00 00 80 41 00 00 88 41  00 00 90 41 00 00 98 41  00 00 00 00 00 00 00 00  00 00 00 00
                  ^ 
    Segmentation fault (core dump)
    
[/code]

alignment についてですが、最大長が `__BIGGEST_ALIGNMENT__` で定義されています。

[code]

    $ gcc -dM -E - < /dev/null | grep __BIGGEST_ALIGNMENT__
    #define __BIGGEST_ALIGNMENT__ 16
    
    $ gcc -mavx2 -dM -E - < /dev/null | grep __BIGGEST_ALIGNMENT__
    #define __BIGGEST_ALIGNMENT__ 32
    
[/code]

AVX2 では 32 bytes(256 bits)なので、gcc の `__attribute__` キーワードによる aligned
属性の指定ができます。

[code]

    float a[64] __attribute__((aligned(32)));
                   :
    
[/code]


CC BY-SA 4.0
2021-01-05T16:04:52.787
73103
2021-01-05T16:04:52.787
73097
answer


回答ではないですが面白かったので参考情報です。

Visual
C++の場合、組み込み関数の意味を踏まえてコード生成してくれます。`prod_avx`のループですが命令の順序が並び変わっていました。あとレジスタをうまくやりくりしてソースコードよりもレジスタ使用数を減らしてもいました。更にコンパイルオプション
`/arch:AVX2` により明示的にAVX2命令の使用を許可したところ勝手にFMA命令まで使って次のようなコードになっていました。

[code]

    ai6 = _mm256_broadcast_ss(a + 8 * i + 6);
    ai7 = _mm256_broadcast_ss(a + 8 * i + 7);
    ai6 = _mm256_mul_ps(ai6, b6n);
    ai4 = _mm256_broadcast_ss(a + 8 * i + 4);
    ai6 = _mm256_fmadd_ps(ai7, b7n, ai6);
    ai5 = _mm256_broadcast_ss(a + 8 * i + 5);
    ai4 = _mm256_mul_ps(ai4, b4n);
    ai2 = _mm256_broadcast_ss(a + 8 * i + 2);
    ai4 = _mm256_fmadd_ps(ai5, b5n, ai4);
    ai3 = _mm256_broadcast_ss(a + 8 * i + 3);
    ai4 = _mm256_add_ps(ai4, ai6);
    ai2 = _mm256_mul_ps(ai2, b2n);
    ai0 = _mm256_broadcast_ss(a + 8 * i);
    ai2 = _mm256_fmadd_ps(ai3, b3n, ai2);
    ai1 = _mm256_broadcast_ss(a + 8 * i + 1);
    ai0 = _mm256_mul_ps(ai0, b0n);
    ai0 = _mm256_fmadd_ps(ai1, b1n, ai0);
    ai0 = _mm256_add_ps(ai0, ai2);
    ai0 = _mm256_add_ps(ai0, ai4);
    
[/code]

質問のソースコードでは、演算順序の都合で全てをFMA命令にすることはできません。

* * *

そんなわけで明示的にFMA命令を使うとこんな感じでしょうか。

[code]

    __m256 ai0, ai1;
    ai0 = _mm256_broadcast_ss(a + 8 * i);
    ai0 = _mm256_mul_ps(ai0, b0n);
    ai1 = _mm256_broadcast_ss(a + 8 * i + 1);
    ai0 = _mm256_fmadd_ps(ai1, b1n, ai0);
    ai1 = _mm256_broadcast_ss(a + 8 * i + 2);
    ai0 = _mm256_fmadd_ps(ai1, b2n, ai0);
    ai1 = _mm256_broadcast_ss(a + 8 * i + 3);
    ai0 = _mm256_fmadd_ps(ai1, b3n, ai0);
    ai1 = _mm256_broadcast_ss(a + 8 * i + 4);
    ai0 = _mm256_fmadd_ps(ai1, b4n, ai0);
    ai1 = _mm256_broadcast_ss(a + 8 * i + 5);
    ai0 = _mm256_fmadd_ps(ai1, b5n, ai0);
    ai1 = _mm256_broadcast_ss(a + 8 * i + 6);
    ai0 = _mm256_fmadd_ps(ai1, b6n, ai0);
    ai1 = _mm256_broadcast_ss(a + 8 * i + 7);
    ai0 = _mm256_fmadd_ps(ai1, b7n, ai0);
    
[/code]

もしくはロードを１回にして、レジスタ上でブロードキャストするとこんな感じでしょうか。

[code]

    __m256 ai0, ai1, ai2, ai3;
    ai3 = _mm256_load_ps(a + 8 * i);                // ai3 = a7, a6, a5, a4, a3, a2, a1, a0
    ai2 = _mm256_permute_ps(ai3, 0x00);             // ai2 = a4, a4, a4, a4, a0, a0, a0, a0
    ai1 = _mm256_permute2f128_ps(ai2, ai2, 0x00);   // ai1 = a0, a0, a0, a0, a0, a0, a0, a0
    ai0 = _mm256_mul_ps(ai1, b0n);                  // ai0 = ai1 * b0n
    ai2 = _mm256_permute2f128_ps(ai2, ai2, 0x11);   // ai2 = a4, a4, a4, a4, a4, a4, a4, a4
    ai0 = _mm256_fmadd_ps(ai2, b4n, ai0);           // ai0 = ai2 * b4n + ai0
    ai2 = _mm256_permute_ps(ai3, 0x55);             // ai2 = a5, a5, a5, a5, a1, a1, a1, a1
    ai1 = _mm256_permute2f128_ps(ai2, ai2, 0x00);   // ai1 = a1, a1, a1, a1, a1, a1, a1, a1
    ai0 = _mm256_fmadd_ps(ai1, b1n, ai0);           // ai0 = ai1 * b1n + ai0
    ai2 = _mm256_permute2f128_ps(ai2, ai2, 0x11);   // ai2 = a5, a5, a5, a5, a5, a5, a5, a5
    ai0 = _mm256_fmadd_ps(ai2, b5n, ai0);           // ai0 = ai2 * b5n + ai0
    ai2 = _mm256_permute_ps(ai3, 0xAA);             // ai2 = a6, a6, a6, a6, a2, a2, a2, a2
    ai1 = _mm256_permute2f128_ps(ai2, ai2, 0x00);   // ai1 = a2, a2, a2, a2, a2, a2, a2, a2
    ai0 = _mm256_fmadd_ps(ai1, b2n, ai0);           // ai0 = ai1 * b2n + ai0
    ai2 = _mm256_permute2f128_ps(ai2, ai2, 0x11);   // ai2 = a6, a6, a6, a6, a6, a6, a6, a6
    ai0 = _mm256_fmadd_ps(ai2, b6n, ai0);           // ai0 = ai2 * b6n + ai0
    ai2 = _mm256_permute_ps(ai3, 0xFF);             // ai2 = a7, a7, a7, a7, a3, a3, a3, a3
    ai1 = _mm256_permute2f128_ps(ai2, ai2, 0x00);   // ai1 = a3, a3, a3, a3, a3, a3, a3, a3
    ai0 = _mm256_fmadd_ps(ai1, b3n, ai0);           // ai0 = ai1 * b3n + ai0
    ai2 = _mm256_permute2f128_ps(ai2, ai2, 0x11);   // ai2 = a7, a7, a7, a7, a7, a7, a7, a7
    ai0 = _mm256_fmadd_ps(ai2, b7n, ai0);           // ai0 = ai2 * b7n + ai0
    
[/code]


CC BY-SA 4.0
2021-01-09T01:00:06.340
73162
2021-01-09T01:00:06.340
4236
73097
answer