裏目小僧の部屋でLazarusの話

Lazarusでx64 SIMD命令を使ってみる

最初に、自分のPCのCPUがどの命令をサポートしてるかを確認します

CPUIDはインラインアセンブラでないと確認出来ないので 呼び出す部分を関数にします

//帰り値用のレコード型
type TCPUIDrec = record  EAX, EBX, ECX, EDX: uint32; end;
//レコード型の帰り値はポインタとなり隠れた引数となります
//つまり procedure getCPUID(var RCX:TCPUIDrec;sEAX,sECX:UINT32); と呼ばれたのと処理は同じです
//CPUIDを呼ぶのにRCX が壊れ CPUIDでまたRCXが壊れるのでpush popにxchgを使っています
{$ASMMODE intel}     
function getCPUID(sEAX,sECX:UINT32): TCPUIDrec; assembler;
  asm //RCX:帰値   RDX  R8d           RCX
  push  RCX    //RCXを保存する
  MOV     EAX,RDX
  MOV     ECX,R8d
  CPUID
  xchg  RCX,[RSP] //RCXを回復すると同時に今のRCXを保存
  MOV   DWORD ptr  0[ RCX],EAX //
  MOV   DWORD ptr  4[ RCX],EBX //
  pop   RAX
  MOV   DWORD ptr  8[ RCX],EAX //  ECX
  MOV   DWORD ptr 12[ RCX],EDX //
end;
type
  TSSEMODE = (sseSSE2, sseSSE3, sseFMA, sseSSE41, sseSSE42, sseAVX, sseAVX2, sseAVX512F, sseAVX512DQ);
  TSSEMODEset = set of TSSEMODE;
const SSEname: array [low(TSSEMODE)..High(TSSEMODE)] of string =
    ('SSE2', 'SSE3', 'FMA', 'SSE41', 'SSE42', 'AVX', 'AVX2', 'AVX512F', 'AVX512DQ');
                                       

function canSIMDused(): TSSEMODEset;
var iSZ, iDX, iCX, iBX: DWORD;
begin
   with  getCPUID(0,0) do iSZ:=EAX;
   with  getCPUID(1,0) do begin iDX:=EAX; iCX:=ECX; end;
   with  getCPUID(7,0) do iBX:=EBX;
//  asm
//           MOV     EAX,0
//           CPUID
//           MOV     iSZ,EAX
//           MOV     EAX,1
//           CPUID
//           MOV     iDX,EDX
//           MOV     iCX,ECX
//           MOV     EAX,7
//           MOV     ECX,0
//           CPUID
//           MOV     iBX,EBX
//  end ['EAX', 'ECX', 'EDX', 'EBX'];

  if iSZ < 7 then iBX := 0;
  Result := [];
  if (iDX and (1 shl 26)) <> 0 then Include(Result, sseSSE2);
  if (iCX and (1 shl 0)) <> 0 then Include(Result, sseSSE3);
  if (iCX and (1 shl 12)) <> 0 then Include(Result, sseFMA);
  if (iCX and (1 shl 19)) <> 0 then Include(Result, sseSSE41);
  if (iCX and (1 shl 20)) <> 0 then Include(Result, sseSSE42);
  if (iCX and (1 shl 28)) <> 0 then Include(Result, sseAVX);
  if (iBX and (1 shl 5)) <> 0 then Include(Result, sseAVX2);
  if (iBX and (1 shl 16)) <> 0 then Include(Result, sseAVX512F);
  if (iBX and (1 shl 17)) <> 0 then Include(Result, sseAVX512DQ);
end;     
  • 私はあまり最新のPCを持っていないので AVX2あたりまでのサポート確認 AVX512は無い事の確認の為に追加

確認方法

  • メモでも張り付けるか コンソールモードで
var d: TSSEMODEset;
 i: TSSEMODE;
 s: string;
begin
 s := '';
 d := GetSSE();
 for i := low(SSEname) to High(SSEname) do
   if i in d then s := s + ' ' + SSEname[i];

で s を表示します。

並列で D[i]:=D[i]+c*S[i] を求める関数 AVXcoffSumShort.inc(31)

  • SSE命令はsingle(32bit浮動小数点)を4個並列に計算出来ます
  • SSE命令はWin64なら必ず使えます。
  • AVXのymmは256bitあり32bitのshortなら8個を一度に求められます。
  • さらにVの付いた命令はSSEと違ってアライメントがそろっていなくても例外が出なく実行されます。
  • もちろんアライメントミスで倍時間はかかるようだけど、それでも並列で処理されるので十分高速になります。
  • ループの処理のオーバーヘッドがあるので 並列処理命令を4つ並べて 32個を1ループで処理します
  • ymmレジスタは16個ありますからもっと増やしても良いのですが・・・
  • canUSEAVX というフラグがfalseだと SSE命令で処理します 上のチェックでAVX+FMAならtrueにして下さい
呼び出す時は arD arSには配列の最初の値(参照なのでアドレス)を渡し coffには数値を渡します
var Dst,Src: array [0..1023] of single;
  coff:single;
//というふうな感じ
 AVXcoffSum(Dst[0], Src[0], coff, length(Dst),true);
  useAVX: boolean はAVX/FMA命令を使う時にtrueに使えない時はfalseにします

引数をポインタにすればいいんだけど、それにはポインタ型を定義しないといけないので 面倒。オープン配列とかにすると 動的配列型や配列型を渡すのにD5だと面倒だったので、このスタイルに。

  • {$ASMMODE intel} はIntelのアセンブラの文法で表現するという指令です。FPC/GCCのLinux系とは転送方向が逆になるので そちらに慣れている方はご注意下さい。
  • SSE命令だと掛け算と加算が必要なのに FMAが有効だと VFMADD231PS と1つの命令で掛算と加算が出来てしましますね
  • AVXcoffSumは引数の並びで決まるレジスタを利用してるので 関数内関数やメソッドにすると隠れた引数が1つ増える為に利用するレジスタを変更する必要があります。

FMA 命令の結果とSSE命令では結果が違う

  • 同じ計算をしてる上の例でもFMA(AVX)を使うかどうかで結果が違います。
    • 当然 SSEによる並列計算とfor文でした結果は同じです。
  • 理由はFMAは掛算と加算を一度にしてるのに対してSSEは掛け算の結果を丸めてから加算するしかないからです
  • 最下位ビットの誤差 単精度で1/2^24 でも6E-8程あるので、結構大きく見えてしまいます

AVXcoffSumShort.inc(31)

ファイルの文字コードはこのページでの表示の為にShiftJisになっています。エディタで変換するか、下の表示コードをIDEにコピーぺしてお使い下さい。

// for i:=0 to  Count-1 do  arD[i]:=arD[i] + coff*arS[i];
{$ASMMODE intel}
{$IFDEF WIN64}
procedure AVXcoffSumShort(pD, pS: pSingle; const coff: single; Count: integer;useAVX: boolean );{$IFDEF FPC} assembler;{$ENDif}
asm //                        RCX  RDX           xmm2          R9
         cmp  useAVX,0
        je  @NoAVX
        VBROADCASTSS ymm2,XMM2  //  xmm2がcoff
        CMP         Count,8
        JMP         @skp1                //while Count >= 32 do
@whilelp24:
        VMOVUPS     YMM0,       [RCX]       //AVX
        VFMADD231PS YMM0, YMM2, [RDX] //FMA YMM0+=YMM5*[RDX]
        VMOVUPS                 [RCX],ymm0
        VMOVUPS     YMM0,       [RCX+32]       //AVX
        VFMADD231PS YMM0, YMM2, [RDX+32] //FMA YMM0+=YMM5*[RDX]
        VMOVUPS                 [RCX+32],    ymm0

        ADD         RCX,40h    //Inc(pD, 8*2)
        ADD         RDX,40h    //Inc(pS, 8*2)
        SUB         Count,16    //Dec(Count, 8*2);
@whilelp8:
        VMOVUPS     YMM0, [RCX]       //AVX
        VFMADD231PS YMM0, YMM2, [RDX] //FMA YMM0+=YMM5*[RDX]
        VMOVUPS     [RCX],    ymm0

        ADD         RCX,20h    //Inc(pD, 8)
        ADD         RDX,20h    //Inc(pS, 8)
        SUB         Count,8    //Dec(Count, 8);
@skp1:
        CMP         Count,24
        JGE         @whilelp24
        CMP         Count,8
        JGE         @whilelp8
        VZEROUPPER        //SSE/AVX切り替えペナルティをなくすために
 @NoAVX:  // 以下 SSE命令
        SHUFPS      XMM2,XMM2,0       //XMM2の最下位で4個埋める
        JMP         @skp2                //while Count >= 32 do
@whilelp12:
        MOVUPS      XMM0, [RDX]       //4word一度に処理する
        MOVUPS      XMM1, [RCX]       //4word一度に処理する
        MULPS       XMM0, XMM2        //coff*pS^
        ADDPS       XMM1, XMM0
        MOVUPS            [RCX],XMM1
        MOVUPS      XMM0, [RDX+16]       //4word一度に処理する
        MOVUPS      XMM1, [RCX+16]       //4word一度に処理する
        MULPS       XMM0, XMM2        //coff*pS^
        ADDPS       XMM1, XMM0
        MOVUPS            [RCX+16],XMM1
        ADD         RCX,32    //Inc(pD, 4*2)
        ADD         RDX,32    //Inc(pS, 4*2)
        SUB         Count,8   //Dec(Count, 4*2);
@whilelp4:
        MOVUPS      XMM0, [RDX]       //4word一度に処理する
        MOVUPS      XMM1, [RCX]       //4word一度に処理する
        MULPS       XMM0, XMM2        //coff*pS^
        ADDPS       XMM1, XMM0
        MOVUPS            [RCX],XMM1
        ADD         RCX,16    //Inc(pD, 4)
        ADD         RDX,16    //Inc(pS, 4)
        SUB         Count,4   //Dec(Count, 4);
@skp2:
        CMP         Count,12
        JGE         @whilelp12
        CMP         Count,4
        JGE         @whilelp4
//残りは0〜3word
          CMP         Count,0
          JL         @skp3

        MOVUPS      XMM0, [RDX]       //4word一度に処理する上位はゴミ
        MOVUPS      XMM1, [RCX]       //
        MULPS       XMM0, XMM2        //coff*pS^
        ADDPS       XMM1, XMM0

        MOVSS            [RCX],XMM1   //最下位保存
        CMP         Count,1
        JL         @skp3
         SHUFPS      XMM1,XMM1,39h   //  4つのデータを右回転=00 11 10 01b
        MOVSS            [RCX+4],XMM1
        CMP         Count,2
        JL         @skp3
        SHUFPS      XMM1,XMM1,39h   //  4つのデータを右回転=00 11 10 01b
        MOVSS            [RCX+8],XMM1
@skp3:
end;
{$ELSE}
procedure AVXcoffSumShort(arD, arS: pSingle; const coff: single; Count: integer;useAVX: boolean );{$IFDEF FPC} assembler;{$ENDif}
asm //                         EAX  RDX                 [$c+ebp]     ecx            [8+ebp]
       MOVSS XMM2,coff;
       cmp  useAVX,0
        je  @NoAVX
        VBROADCASTSS ymm2,XMM2  //  xmm2がcoff
        CMP         Count,8
        JMP         @skp1                //while Count >= 32 do
@whilelp24:
        VMOVUPS     YMM0,       [EAX]       //AVX
        VFMADD231PS YMM0, YMM2, [EDX] //FMA YMM0+=YMM5*[EDX]
        VMOVUPS                 [EAX],ymm0
        VMOVUPS     YMM0,       [EAX+32]       //AVX
        VFMADD231PS YMM0, YMM2, [EDX+32] //FMA YMM0+=YMM5*[EDX]
        VMOVUPS                 [EAX+32],    ymm0

        ADD         EAX,40h    //Inc(pD, 8*2)
        ADD         EDX,40h    //Inc(pS, 8*2)
        SUB         Count,16    //Dec(Count, 8*2);
@whilelp8:
        VMOVUPS     YMM0, [EAX]       //AVX
        VFMADD231PS YMM0, YMM2, [EDX] //FMA YMM0+=YMM5*[EDX]
        VMOVUPS     [EAX],    ymm0

        ADD         EAX,20h    //Inc(pD, 8)
        ADD         EDX,20h    //Inc(pS, 8)
        SUB         Count,8    //Dec(Count, 8);
@skp1:
        CMP         Count,24
        JGE         @whilelp24
        CMP         Count,8
        JGE         @whilelp8
        VZEROUPPER       //SSE/AVX切り替えペナルティをなくすために
@NoAVX:  // 以下 SSE命令
        SHUFPS      XMM2,XMM2,0       //XMM2の最下位で4個埋める
        MOV         EDX,arS
        MOV         EAX,arD
        JMP         @skp2                //while Count >= 32 do
@whilelp12:
        MOVUPS      XMM0, [EDX]       //4word一度に処理する
        MOVUPS      XMM1, [EAX]       //4word一度に処理する
        MULPS       XMM0, XMM2        //coff*pS^
        ADDPS       XMM1, XMM0
        MOVUPS            [EAX],XMM1
        MOVUPS      XMM0, [EDX+16]       //4word一度に処理する
        MOVUPS      XMM1, [EAX+16]       //4word一度に処理する
        MULPS       XMM0, XMM2        //coff*pS^
        ADDPS       XMM1, XMM0
        MOVUPS            [EAX+16],XMM1
        ADD         EAX,32    //Inc(pD, 4*2)
        ADD         EDX,32    //Inc(pS, 4*2)
        SUB         Count,8   //Dec(Count, 4*2);
@whilelp4:
        MOVUPS      XMM0, [EDX]       //4word一度に処理する
        MOVUPS      XMM1, [EAX]       //4word一度に処理する
        MULPS       XMM0, XMM2        //coff*pS^
        ADDPS       XMM1, XMM0
        MOVUPS            [EAX],XMM1
        ADD         EAX,16    //Inc(pD, 4)
        ADD         EDX,16    //Inc(pS, 4)
        SUB         Count,4   //Dec(Count, 4);
@skp2:
        CMP         Count,12
        JGE         @whilelp12
        CMP         Count,4
        JGE         @whilelp4
//残りは0〜3word
          CMP         Count,0
          JL         @skp3
        MOVUPS      XMM0, [EDX]       //4word一度に処理する上位はゴミ
        MOVUPS      XMM1, [EAX]       //
        MULPS       XMM0, XMM2        //coff*pS^
        ADDPS       XMM1, XMM0

        MOVSS            [EAX],XMM1   //最下位保存
        CMP         Count,1
        JL         @skp3
         SHUFPS      XMM1,XMM1,39h   //  4つのデータを右回転=00 11 10 01b
        MOVSS            [EAX+4],XMM1
        CMP         Count,2
        JL         @skp3
        SHUFPS      XMM1,XMM1,39h   //  4つのデータを右回転=00 11 10 01b
        MOVSS            [EAX+8],XMM1
@skp3:
end;
{$ENDIF}
//テスト長の同値を返す関数(useAVXがtrueなら近似値になる)
procedure testAVXcoffSumShort( arD, arS: psingle; const coff: single; Count: integer;useAVX: boolean );
var i:Integer;
begin
for i:=0 to  Count-1 do
  arD[i]:=arD[i] + coff*arS[i];
end;

プライバシーポリシー本文は日本語以外に翻訳禁止 click寄付