SIMDプログラミングによる行列積について


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]




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]

