3

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

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

以下がコードです。

#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);
    
}
2
  • 1
    少なくともメモリアラインメントの問題があるように思われます。コンパイラ、プラットフォームは何ですか?
    – OOPer
    Commented 2021年1月5日 12:51
  • OS:ubuntu 18.04、カーネル:4.4.0-143-generic、x64、コンパイラはgcc version 5.4.0 です。
    – puzzlesimd
    Commented 2021年1月5日 13:23

3 件の回答 3

2

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 を指定しています)。

$ 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)

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

$ 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

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

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

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

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

prod_avx内のcdの宣言:

    float *cd = (float *)aligned_alloc(32, 64 * sizeof(float));

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

    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));

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

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

0

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

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

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);

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


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

__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);

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

__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

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。