Upgrade to Pro — share decks privately, control downloads, hide ads and more …

プログラムを高速化する話

 プログラムを高速化する話

プログラムを高速化するためのテクニックをまとめました。GPU編はこちら https://speakerdeck.com/primenumber/puroguramuwogao-su-hua-suruhua-ii-gpgpubian

prime number

July 10, 2023
Tweet

More Decks by prime number

Other Decks in Programming

Transcript

  1. 2 自己紹介 第 37 代 KMC 会計( 2014/02 〜 2015/01

    ) root ( 2013/11 〜) 理学部数学系志望 ( もうすぐ系登録試験 ) 競技プログラミング練習会 GR-SAKURA で遊ぼう コーディング大海 KMC 学習発表会
  2. 3 自己紹介 過去の KMC での活動 難解プログラミング言語勉強会 オセロの AI を作ろう 2048AI

    コンテスト RICOH and Java Developer Challange Plus 2013/2014 KMC 流しそうめん KMC 静止そうめん KMC 鍋 prime/ 免許合宿
  3. 11 最適化について 「細かい効率のことは忘れて、時間の 97% について考え よう。時期尚早な最適化は諸悪の根源だ。それでも残り 3% についても機会を逃すべきではない」 - Donald

    E. Knuth 「プログラム最適化の第一法則 : 最適化するな。 プログラム最適化の第二法則 ( 上級者限定 ): まだするな。 」 - Michael A. Jackson
  4. 13 最適化の対象 主に Intel の Haswell マイクロアーキテクチャ以降を対象 多くのテクニックは他のプロセッサにも応用できます ベース マイクロアーキテクチャ

    プロセスルール 登場年 Nehalem Nehalem 45nm 2008 〃 Westmere 32nm 2010 Sandy Bridge Sandy Bridge 32nm 2011 〃 Ivy Bridge 22nm 2012 Haswell Haswell 22nm 2013 〃 Broadwell 14nm 2014
  5. 15 Intel C/C++ Compiler Intel 謹製の C/C++ コンパイラ 最適化能力が強いと言われている 学生なら非営利目的に限り無料で入手可能

    Intel Parallel Studio XE についてくる 最新版は gcc 4.9 互換なので C++ の最新規格のサポートも ぼちぼち
  6. 18 最適化手法を学ぶには 「インテル ® 64 アーキテクチャー および IA-32 アーキ テクチャー

    最適化リフ レンス マニュアル」 [0][1] を読め (完) PDF ファイルをインテルの公式サイトからダウンロード できます 日本語訳もあるけどちょっと古い 700 ページ以上の超大作
  7. 23 メモリアクセスの遅延 転送速度も足りない メモリ側は DDR4-4266 でも 34.1GB/s CPU 側は最大で読み込み 64Bytes/cycle

    、書き込み 32Bytes/cycle なので、 CPU の動作周波数が 3GHz なら、 読み込み 192GB/s 、書き込み 32Bytes/cycle 線形にアクセスする場合でも速度が足りない!
  8. 26 局所的でないメモリアクセスを避ける 局所的でないメモリアクセスをすると、アクセスする データがキャッシュに乗っている確率が低くなる 例 ) 行列積 for (int i

    = 0; i < ROWS; ++i) for (int j = 0; j < COLS; ++j) for (int k = 0; k < LEN; ++k) C[i][j] += A[i][k] * B[k][j]; 二次元配列 B に ROWS 要素飛びでアクセスしている
  9. 27 局所的でないメモリアクセスを避ける 局所的でないメモリアクセスをすると、アクセスする データがキャッシュに乗っている確率が低くなる 例 ) 行列積 for (int i

    = 0; i < ROWS; ++i) for (int k = 0; k < LEN; ++k) for (int j = 0; j < COLS; ++j) C[i][j] += A[i][k] * B[k][j]; 全ての配列に順番にアクセスするようになった 入れ替えた
  10. 28 データ構造を SoA に 大量のデータを順番に処理するとき、 AoS(Array of Structs; 構造体の配列 )

    よりも、 SoA(Struct of Arrays; 配列の構造体 ) の方が高速に動作 する可能性がある
  11. 29 データ構造を SoA に 例 struct data { int a,

    b, c; double x, y, z; } d_ary[SIZE]; // AoS int a[SIZE], b[SIZE], c[SIZE]; double x[SIZE], y[SIZE], z[SIZE]; //SoA
  12. 30 データ構造を SoA に a0 b0 c0 x0 y0 z0

    a1 b1 c1 x1 y1 z1 a2 b2 c2 … SoAだと順番にaにアクセスすると6要素ごとになる AoSだと順番にaにアクセスすると 連続した領域にアクセスできる SIMD命令を使いやすくなる a0 a1 a2 … b0 b1 b2 … c0 c1 c2 … x0 x1 x2 … y0 y1 y2 … z0 z1 z2 …
  13. 32 ストリップマイニング 以下のコードを考える for (int i = 0; i <

    SIZE; ++i) { hoge(A[i]); } for (int i = 0; i < SIZE; ++i) { fuga(A[i]); } 配列 A が十分長いとき、最初のループが終わった時点で、 A の先頭はキャッシュから排出されている
  14. 33 ストリップマイニング したがって、 A は 2 回メインメモリから読み込まれるこ とになり、効率が悪い 一方で、 hoge

    や fuga は SIMD 命令で並列に処理できる ものとする ここで、ループをキャッシュに乗るサイズに分割すると、 SIMD 命令を使いつつ、メインメモリへのアクセスを 1 回 に減らせる
  15. 34 ストリップマイニング for (int i = 0; i < SIZE;

    i += strip_size) { for (int j = i; j < min(SIZE, i+strip_size); ++j) { hoge(A[j]); } for (int j = i; j < min(SIZE, i+strip_size); ++j) { fuga(A[j]); } }
  16. 35 ブロック化 次のコードを考える for (int i = 0; i <

    SIZE; ++i) { for (int j = 0; j < SIZE; ++j) { A[i][j] += B[j][i]; } } B には飛び飛びのアクセスをしているので SIZE が大きい と毎回キャッシュミスが発生して効率が悪い
  17. 38 ブロック化 for (int i = 0; i < SIZE;

    i+=block_size) { for (int j = 0; j < SIZE; j+=block_size) { for (int ii = i; ii < i+block_size; ++ii) { for(int jj = j; jj < j+block_size; ++jj) { A[ii][jj] += B[jj][ii]; } } } }
  18. 40 基礎知識 コンピュータ内部ではデータは 2 進数で管理されている 2 進数リテラルは C++ では 0b

    のあとに続けて 0/1 を書く ことで記述できる 1 0 1 1 1 0 0 1 1バイト = 8ビット 0b10111001
  19. 41 基礎知識 符号なし整数は単純に 2 進数で表される 1 バイトなら 0 〜 255=2

    -1 ⁸ 4 バイトなら 0 〜 4294967295=2³²-1 n バイトなら 0 〜 2^8n-1 0b10111001=185
  20. 42 基礎知識 符号付き整数は一番上の桁が符号を表す ( 負の数なら 1) n ビット符号付き整数は、正の数なら符号なしと同じ 負の数なら、 mod

    2^n で同じになる正の数の 2 進数表記 0b10111001 を 1 バイト符号付き整数だとみなすと、 0b10111001+0b01000111=2⁸ なので、 0b10111001=-0b01000111=-71
  21. 43 基礎知識 符号付き整数は一番上の桁が符号を表す ( 負の数なら 1) n ビット符号付き整数は、正の数なら符号なしと同じ 負の数なら、 mod

    2^n で同じになる正の数の 2 進数表記 1 バイトなら -128 〜 127 4 バイトなら -2147483648 〜 2147483647 n バイトなら -2^(n-1) 〜 2^(n-1)-1
  22. 44 基礎知識 実数は 1.x×2^i の形の有理数に丸めて表現する 丸める精度に応じてバリエーションがある 単精度浮動小数型 4 バイト、 2

    進数で 23 桁の精度、 ±10^±38 程度まで表現 できる 倍精度浮動小数型 8 バイト、 2 進数で 53 桁の精度、 ±10^±308 程度まで表 現できる
  23. 45 ビット演算 2 進数の 0/1 の列を操作するような演算の総称 ビット論理和 ビット論理積 ビット排他的論理和 ビット否定

    ビットシフト 加減乗算などがビット列を操作するのに使われることも 最近はビット列操作用の命令も追加されている
  24. 48 基本的なビット演算 ビット論理和 OR (C 言語 : |) 片方でも 1

    ならば 1 、そうでないなら 0 0 0 1 0 1 1 0 0 1 0 1 1 1 0 0 1 1 0 1 1 1 1 0 1 A B A OR B
  25. 49 基本的なビット演算 ビット論理積 AND (C 言語 : &) 両方とも 1

    ならば 1 、そうでないなら 0 0 0 1 0 1 1 0 0 1 0 1 1 1 0 0 1 0 0 1 0 1 0 0 0 A B A AND B
  26. 50 基本的なビット演算 ビット排他的論理和 XOR (C 言語 : ^) 片方だけ 1

    ならば 1 、そうでないなら 0 0 0 1 0 1 1 0 0 1 0 1 1 1 0 0 1 1 0 0 1 0 1 0 1 A B A XOR B
  27. 52 基本的なビット演算 ビットシフト (C 言語 : 左シフト <<, 右シフト >>)

    0/1 列を右や左にシフトする 左シフトが上の桁の方にシフトする 1 1 1 0 0 1 0 0 1 0 1 1 1 0 0 1 A A << 2 シフトしたときに詰める数字によっていくつかバリ エーションが存在する
  28. 53 C 言語でのビット演算 // ビット演算は整数型に対してのみ使える A = B & C;

    A &= B; A = B << 4; // 下位ビットには 0 が詰められる A <<= 4; A = ~B; // 符号なし 64 ビットリテラルを扱うとき A = UINT64_C(0xCCCCCCCCCCCCCCCC);
  29. 55 特定のビットを操作する UINT64_C(1) << index で index ビット目のみ立った数を 表せる //A

    の index ビット目を操作する A |= UINT64_C(1) << index; // ビットを立てる A &= ~(UINT64_C(1) << index); // ビットを下ろす A ^= UINT64_C(1) << index; // ビットを反転する result = (A >> index) & 1; // ビットの値を取得する
  30. 56 マスク // 奇数ビット目をクリア // 0x5 = 0b0101 A &=

    UINT64_C(0x5555555555555555); // 偶数ビットなら 0xAAA... = 0b1010... // 2 ビットごとに交互にクリアするなら 0xCCC...=0b11001100...
  31. 59 立っているビットの数を数える (popcount) 1 0 1 1 1 0 0

    1 1 1 1 0 1 1 1 0 A A&0xAA (A&0xAA) >> 1 0 1 0 1 A&0x55 0 1 1 0 0 1 0 1 ((A&0xAA) >> 1) + (A&0x55) 2ビットごとの立っているビットの数の和
  32. 60 立っているビットの数を数える (popcount) 0 1 1 0 0 1 0

    1 0 1 0 1 0 1 0 1 A' A'&0xCC (A'&0xCC) >> 2 1 0 0 1 A'&0x33 0 0 1 1 0 0 1 0 ((A'&0xCC) >> 2) + (A'&0x33) 4ビットごとの立っているビットの数の和
  33. 61 立っているビットの数を数える (popcount) 0 0 1 1 0 0 1

    0 0 0 1 1 0 0 1 1 A'' A''&0xF0 (A''&0xF0) >> 4 0 0 1 0 A''&0x0F 0 0 0 0 0 1 0 1 ((A''&0xF0) >> 4) + (A''&0x0F) 8ビット全体の立っているビットの数の和
  34. 64 複数のビットから成るデータの配列の  ハミング距離を求める 例えば、 2 ビットから成るデータの配列のハミング距離 を求めるとき、 for (int i

    = 0; i < ARY_SIZE; ++i) { uint8_t C = A[i] ^ B[i]; //8 ビットの場合 C = ((C & 0xAA) >> 1) | (C & 0x55); result += popcount(C); }
  35. 65 立っている一番下のビットを求める B = A & -A; 1 0 1

    1 1 0 0 0 A 0 1 0 0 1 0 0 0 -A 0 0 0 0 1 0 0 0 A & -A -A = ~A + 1であることを利用
  36. 66 立っている一番下のビットを中心に操作 A & (A – 1); // 立っている一番下のビットをクリア A

    ^ -A; // 立っている一番下のビットより上の桁を 1 に A | -A; // さらに立っている一番下のビットも 1 に // 立っている一番下のビットより下の桁を 1 に A ^ (A – 1)
  37. 67 立っているビット列を走査する // i &= i-1 で i の立っている一番下のビットをクリア for

    (uint64_t i = bits; i != 0; i &= i-1) { uint64_t rmb = i & -i; // 何らかの処理 } 立っているビットの数が少ない場合には、この方法でも 立っているビットの数を高速に数えられる
  38. 68 立っている一番上のビットを求める 二分探索で求める 0 0 1 1 1 0 0

    1 A 0 0 1 1 A&0xF0 != 0 立っている一番上のビットは上4桁にある
  39. 71 立っている一番上のビットを求める // 8 ビットの場合 A = (A & 0xF0)

    ? (A & 0xF0) : A; A = (A & 0xCC) ? (A & 0xCC) : A; A = (A & 0xAA) ? (A & 0xAA) : A;
  40. 74 ビット列の並びを反転する 分割統治法を用いる 1 0 1 1 1 0 0

    1 A 1 0 0 1 1 0 1 1 B=(A>>4) | (A<<4) 0 1 1 0 1 1 1 0 ((B&0xCC)>>2) | ((B&0x33)<<2)
  41. 75 ビット列の並びを反転する 分割統治法を用いる 1 0 1 1 1 0 0

    1 A 1 0 0 1 1 0 1 1 B=(A>>4) | (A<<4) 0 1 1 0 1 1 1 0 C=((B&0xCC)>>2) | ((B&0x33)<<2) ((C&0xAA)>>1) | ((C&0x55)<<1) 1 0 0 1 1 1 0 1
  42. 77 ビット列を一部だけスワップする Delta Swap という手法 長さが等しく、重複のないビット列をスワップする * * * A

    B C * * * a b c * * * * この幅をdeltaとする 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 x mask ABCとabcをスワップする deltaとmaskを予め求めておく
  43. 78 ビット列を一部だけスワップする Delta Swap という手法 長さが等しく、重複のないビット列をスワップする * * * A

    B C * * * a b c * * * * 0 0 0 0 0 0 0 0 0 P Q R 0 0 0 0 x b := (x ^ (x >> delta)) & mask P = A^a, Q = B^b, R = C^cとなる b
  44. 79 ビット列を一部だけスワップする Delta Swap という手法 長さが等しく、重複のないビット列をスワップする * * * A

    B C * * * a b c * * * * 0 0 0 P Q R 0 0 0 P Q R 0 0 0 0 x ABCとabcが入れ替わった(他のビットは変化なし) c := b ^ (b << delta) c * * * a b c * * * A B C * * * * c ^ x
  45. 80 ビット列の指定した場所を詰めて並べる 64bit 整数で 8x8 の 2 次元データを表現するとする オセロやチェスなどの盤面の状態、など 0

    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
  46. 83 magic bitboard 実はビットの配置ごとに適切に選んだ整数 magic を掛け てやると 1 列に並ぶ magic

    bitboard として知られている 14 17 21 26 28 42 44 49 53 17 26 44 53 49 42 28 21 14 * magic ) >> some = (
  47. 86 bts, btr, btc, bt 特定のビットを操作する命令 bts unsigned char _bittestandset(

    __int32* a, __int32 b); unsigned char _bittestandset64( __int64* a, __int64 b); 特定のビットを立て、立てる前のそのビットの状態を返す btr: ビットを下ろす btc: ビットを反転する bt: ビットを取得する
  48. 87 blsi, blsmsk, blsr, tzcnt 立っている一番下のビット関連の操作をする命令 blsi: A & -A

    と同じ。立っている一番下のビットを求める blsmsk: A ^ -A と同じ。立っている一番下のビットより 上のビットを 1 にした数を求める blsr: A & (A-1) と同じ。立っている一番下のビットを 0 にした数を求める。 tzcnt: 立っている一番下のビットの桁数を求める
  49. 89 bextr bextr start, len を指定して図のような操作をした結果を返す 0 1 0 0

    1 0 1 1 0 1 0 1 1 1 1 0 start len 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 len
  50. 92 SIMD 命令による最適化 SIMD 命令とは SIMD 命令が使える場合 SIMD 命令の使い方 SIMD

    プログラミングにおけるポイント メモリアラインメント マスクによる条件分岐の除去 その他 SIMD 命令で可能なこと AVX-512 について
  51. 93 SIMD 命令とは SIMD(Single Instruction Multiple Data の略 ) 一つの命令で複数のデータを処理する

    AVX 命令セットでは一度に 256bit(32Bytes) のデータを 一括で処理できる 単精度浮動小数 x8 倍精度浮動小数 x4 1/2/4/8Bytes 整数 x32/16/8/4
  52. 94 SIMD 命令とは A[0] A[1] A[2] A[3] B[0] B[1] B[2]

    B[3] A[0]+B[0] A[1]+B[1] A[2]+B[2] A[3]+B[3] + 256ビット
  53. 97 自動ベクトル化 コンパイラによっては適切なコンパイルオプションを指 定することによって、演算を SIMD 命令に変換してくれる icc なら -xAVX を指定すると

    AVX 命令を使ってくれる しかし、どんな場合も SIMD 命令に変換してくれるわけで はない そのときは、残り 3 通りのいずれかの方法を用いる
  54. 99 SIMD 組み込み関数の利用 #include “immintrin.h” で使えるようになる 単精度浮動小数 x8 型 __m256

    倍精度浮動小数 x4 型 __m256d 整数 x4/8/16/32 型 __m256i __m256d resv = _mm256_add_pd(a, b); のように使う 今回は主に SIMD 組み込み関数を利用していきます
  55. 101 SIMD 組み込み関数の利用例 配列 A の各要素に配列 B の対応する要素を加算した結果 を配列 C

    に格納する。 配列の長さは 4 の倍数、 A,B,C はすべて倍精度浮動小数型 for (int i = 0; i < SIZE/4; ++i) { __m256d va = _mm256_load_pd(A+4*i); __m256d vb = _mm256_load_pd(B+4*i); __m256d res = _mm256_add_pd(va, vb); _mm256_store_pd(C+4*i, res); }
  56. 102 SIMD 組み込み関数の使い方の基本 メモリ上からデータを取ってくるには _mm256_load* メモリ上にデータを書き込むには _mm256_store* _mm256_( 動作 )_(

    型 ) となっていることが多い 型には ps: 単精度浮動小数 x8 pd: 倍精度浮動小数 x4 epi8/16/13/64: 8/16/32/64 バイト整数 x32/16/8/4
  57. 103 Intel Intrinsics Guide https://software.intel.com/sites/landingpage/Intrinsic sGuide/ SIMD 組み込み関数を中心とした各種 CPU 命令の組み込

    み関数のリファレンス 拡張の種類や命令のタイプで絞ったり検索したりできる 便利!!!
  58. 106 メモリアラインメント 動的確保 メモリを動的確保する場合は __alignas は使えない void *_mm_malloc(size_t size, size_t

    align); size バイトの領域を確保。返されるポインタのアドレス は align の倍数になる _mm_malloc で確保したメモリは _mm_free で開放する 他にも posix_memalign/aligned_alloc を使う方法や std::align を使う方法もある
  59. 107 マスクによる条件分岐の除去 __m256i _mm256_cmpeq_epi8(__m256i, __m256i); 各バイトを比較して等しいなら FF 、異なるなら 00 を各

    バイトに入れて返す eq (等しい) /gt (大きい)、 epi8/16/32/64 ( 1/2/4/8 バイト)のバリエーションが存 在する
  60. 110 マスクによる条件分岐の除去 例 ) for (int i = 0; i

    < SIZE; ++i) { if (a[i] > b[i]) { c[i] += a[i]; } } 条件分岐が入っているのでこのままでは並列に足し算で きない
  61. 111 マスクによる条件分岐の除去 例 ) for (int i = 0; i

    < SIZE/32; ++i) { __m256i va = _mm256_load_si256((__m256i*)(A+32*i)); __m256i vb = _mm256_load_si256((__m256i*)(B+32*i)); __m256i mask = _mm256_cmpgt_epi8(va, vb); va = _mm256_and_si256(va, mask); __m256i vc = _mm256_load_si256((__m256i*)(C+32*i)); vc = _mm256_add_epi8(vc, va); _mm256_store_si256((__m256i*)(C+32*i), vc); }
  62. 112 マスクによる条件分岐の除去 例 ) for (int i = 0; i

    < SIZE/32; ++i) { __m256i va = _mm256_load_si256((__m256i*)(A+32*i)); __m256i vb = _mm256_load_si256((__m256i*)(B+32*i)); __m256i mask = _mm256_cmpgt_epi8(va, vb); va = _mm256_and_si256(va, mask); __m256i vc = _mm256_load_si256((__m256i*)(C+32*i)); vc = _mm256_add_epi8(vc, va); _mm256_store_si256((__m256i*)(C+32*i), vc); }
  63. 113 マスクによる条件分岐の除去 例 ) for (int i = 0; i

    < SIZE/32; ++i) { __m256i va = _mm256_load_si256((__m256i*)(A+32*i)); __m256i vb = _mm256_load_si256((__m256i*)(B+32*i)); __m256i mask = _mm256_cmpgt_epi8(va, vb); va = _mm256_and_si256(va, mask); __m256i vc = _mm256_load_si256((__m256i*)(C+32*i)); vc = _mm256_add_epi8(vc, va); _mm256_store_si256((__m256i*)(C+32*i), vc); }
  64. 114 マスクによる条件分岐の除去 例 ) for (int i = 0; i

    < SIZE/32; ++i) { __m256i va = _mm256_load_si256((__m256i*)(A+32*i)); __m256i vb = _mm256_load_si256((__m256i*)(B+32*i)); __m256i mask = _mm256_cmpgt_epi8(va, vb); va = _mm256_and_si256(va, mask); __m256i vc = _mm256_load_si256((__m256i*)(C+32*i)); vc = _mm256_add_epi8(vc, va); _mm256_store_si256((__m256i*)(C+32*i), vc); }
  65. 115 マスクによる条件分岐の除去 例 ) for (int i = 0; i

    < SIZE/32; ++i) { __m256i va = _mm256_load_si256((__m256i*)(A+32*i)); __m256i vb = _mm256_load_si256((__m256i*)(B+32*i)); __m256i mask = _mm256_cmpgt_epi8(va, vb); va = _mm256_and_si256(va, mask); __m256i vc = _mm256_load_si256((__m256i*)(C+32*i)); vc = _mm256_add_epi8(vc, va); _mm256_store_si256((__m256i*)(C+32*i), vc); }
  66. 117 AVX-512 について Broadwell の次の Skylake で Xeon( サーバー向けプロ セッサ

    ) に AVX-512 拡張命令が追加される予定 一度に 512 ビット扱えるだけでなく、ほぼすべての命令 にマスクを掛けることが出来るようになる 各種の便利命令がてんこ盛り Core i3/5/7 等には Skylake の次の Cannonlake で入る?
  67. 120 参考文献等 [0] インテル ® 64 アーキテクチャーおよび IA-32 アーキテク チャー最適化リファレンス・マニュアル

    http://www.intel.co.jp/content/dam/www/public/ijkk/jp/ja/do cuments/developer/248966-024JA.pdf [1] 英語の最新版 http://www.intel.co.jp/content/dam/www/public/us/en/docum ents/manuals/64-ia-32-architectures-optimization-manual.pdf [2] Intel Intrinsics Guide https://software.intel.com/sites/landingpage/Intrin sicsGuide/
  68. 121 参考文献等 [3] Intel® Architecture Instruction Set Extensions Programming Reference

    https://software.intel.com/sites/default/files/managed/0d/53/ 319433-022.pdf [4] Intel® 64 and IA-32 Architectures Software Developer’s Manual http://www.intel.co.jp/content/dam/www/public/us/en/docu ments/manuals/64-ia-32-architectures-software-developer- manual-325462.pdf [5] CPU – Wikipedia http://ja.wikipedia.org/wiki/CPU