画像処理アルゴリズムの高速化2 / CUDA Acceleration Seminar5 20220929

2022年9月29日に開催した「CUDA高速化セミナーvol.5 / 画像処理アルゴリズムの高速化2」の当日資料です。

Transcript

1. Copyright© Fixstars Group CUDA高速化セミナー vol.5 画像処理アルゴリズムの高速化2
2. Copyright© Fixstars Group 本日のAgenda • はじめに • フィックスターズの紹介 • 画像処理アルゴリズムの高速化1 のおさらい • バイラテラルフィルタ高速化 • 転置 • リダクション • まとめ 2
3. Copyright© Fixstars Group はじめに
4. Copyright© Fixstars Group 本講演の位置づけ • CUDAに関連する様々な技術情報を、CUDA高速化セミナーとして発信しています • 今回は、vol.1で解説したデータ転送・カーネルを書く上での基本を踏まえ、CUDA 特有の 計算方法を使った、バイラテラルフィルタ、転置、リダクションの実装例を紹介します • こんな方に向いています ◦ これから CUDA を使った画像処理をしてみたい ◦ CUDA カーネルを高速化したい 4 • vol.1 画像処理アルゴリズムの高速化 • vol.2 CUDAアーキテクチャの進化 • vol.3 ソフトウェア高速化と深層学習 • vol.4 TensorRT化のワークフロー事例紹介
5. Copyright© Fixstars Group 発表者紹介 冨田 明彦 ソリューションカンパニー 執行役員 2008年に入社。金融、医療業界において、 ソフトウェア高速化業務に携わる。その後、 新規事業企画、半導体業界の事業を担当し、 現職。 5 上野 晃司 ソリューション第一事業部 エグゼクティブエンジニア 2016年に入社。スパコンのベンチマーク Graph500を「京」「富岳」向けに最適化し 世界１位を達成。CUDAやOpenCLを使った 画像処理高速化を担当。
6. Copyright© Fixstars Group フィックスターズの ご紹介
7. Copyright© Fixstars Group フィックスターズの強み コンピュータの性能を最大限に引き出す、ソフトウェア高速化のエキスパート集団 ハードウェアの知見 アルゴリズム実装力 各産業・研究分野の知見 7 目的の製品に最適なハードウェアを見抜 き、その性能をフル活用するソフトウェア を開発します。 ハードウェアの特徴と製品要求仕様に合わ せて、アルゴリズムを改良して高速化を実 現します。 開発したい製品に使える技術を見抜き、実 際に動作する実装までトータルにサポート します。
8. Copyright© Fixstars Group サービス概要 お客様専任のエンジニアが直接ヒアリングを行い、高速化を実現するために乗り越えるべき 課題や問題を明確にしていきます。 高速化のワークフロー コンサルティング 先行技術調査 性能評価・ボトルネックの特定 高速化 アルゴリズムの改良・開発 ハードウェアへの最適化 レポート作成 サポート レポートやコードへのQ&A 実製品への組込み支援 8
9. Copyright© Fixstars Group サービス提供分野 9 半導体 自動車 産業機器 生命科学 金融 • NAND型フラッシュメモリ向け ファームウェア開発 • 次世代AIチップの開発環境基盤 • 自動運転の高性能化、実用化 • 次世代パーソナルモビリティの 研究開発 • Smart Factory実現への支援 • マシンビジョンシステムの高速化 • ゲノム解析の高速化 • 医用画像処理の高速化 • AI画像診断システムの研究開発 • デリバティブシステムの高速化 • HFT(アルゴリズムトレード)の高速化
10. Copyright© Fixstars Group サービス領域 様々な領域でソフトウェア高速化サービスを提供しています。大量データの高速処理は、 お客様の製品競争力の源泉となっています。 10 組込み高速化 画像処理・アルゴリズム 開発 分散並列システム開発 GPU向け高速化 FPGAを活用した システム開発 量子コンピューティング AI・深層学習 自動車向け ソフトウェア開発 フラッシュメモリ向け ファームウェア開発
11. Copyright© Fixstars Group 画像処理アルゴリズム開発 高速な画像処理需要に対して、経験豊富なエンジニアが 責任を持って製品開発をご支援します。 お客様の課題 高度な画像処理や深層学習等のアルゴリズム を開発できる人材が社内に限られている 機能要件は満たせそうだが、ターゲット機器 上で性能要件までクリアできるか不安 製品化に結びつくような研究ができていない ご支援内容 深層学習ネットワーク精度の改善 様々な手法を駆使して深層学習ネットワークの精度を改善 論文調査・改善活動 論文調査から最先端の手法の探索 性能向上に向けた改善活動を継続 アルゴリズム調査・改変 課題に合ったアルゴリズム・実装手法を調査 製品実装に向けて適切な改変を実施 11
12. Copyright© Fixstars Group GPU向け高速化 高性能なGPUの本来の性能を十分に引き出し、 ソフトウェアの高速化を実現します。 お客様の課題 GPUで計算してみたが期待した性能が出ない GPU/CPUを組み合わせた全体として最適な設 計がしたい ご支援内容 GPU高速化に関するコンサルティング CPU・GPU混在環境でのシステム設計 アルゴリズムのGPU向け移植 GPUプログラム高速化 継続的な精度向上 原価を維持したまま機能を追加するため、も う少し処理を速くしたい 品質確保のため、精度を上げたく演算量は増 えるが性能は維持したい 12
13. Copyright© Fixstars Group AI・深層学習向け技術支援 AIを使うためのハードウェア選定や、高速な計算を実現する ソフトウェア開発技術で、お客様の製品開発を支援します。 お客様の課題 推論精度を維持したまま計算時間を短縮 したい 組込みデバイス向けにAIモデルを軽量化 したい ご支援内容 AIモデル設計 データの前処理・後処理 推論精度の改善 分散処理による学習高速化 モデル圧縮・推論の高速化 学習計算を高速化して研究開発を効率化 したい 精度と計算時間を両立するAIモデルを 開発したい 13
14. Copyright© Fixstars Group 前回の復習
15. Copyright© Fixstars Group NVIDIA Nsight Compute • カーネルプロファイラをサポート 15
16. Copyright© Fixstars Group ガウシアンフィルタ ⊗ カーネル 16
17. Copyright© Fixstars Group ガウシアンフィルタCUDA化 スレッド割り当て ブロック (0,0) ブロック (0,1) ブロック (0,2) ブロック (0,3) ブロック (0,4) ブロック (1,0) ブロック (1,1) ブロック (1,2) ブロック (1,3) ブロック (1,4) ブロック (2,0) ブロック (2,1) ブロック (2,2) ブロック (2,3) ブロック (2,4) ブロック (3,0) ブロック (3,1) ブロック (3,2) ブロック (3,3) ブロック (3,4) 32 32 • 1スレッドが出力1 ピクセルを担当 • ブロックの最大ス レッド数は1024な ので、１ブロック 32x32(=1024ス レッド)に設定 • 画像全体を覆うよう にブロックを起動す る 17
18. Copyright© Fixstars Group 本日説明するコード • ↓ここにあります • https://github.com/fixstars/CudaOptimizeSample/blob/master/CudaO ptimizeSample/kernel.cu • （前回と同じファイルです） 18
19. Copyright© Fixstars Group バイラテラルフィルタ 高速化
20. Copyright© Fixstars Group バイラテラルフィルタ高速化 • ガウシアンフィルタと同じように、画像と重みを掛け合わせるが、重みが画 素値の差によって変わるようにしたもの 20
21. Copyright© Fixstars Group __global__ void BilateralKernelSimple( const uint8_t *src, uint8_t *dst, int width, int height, int step, float sigma) { int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; if (x < width && y < height) { float coef = 1.0 / sqrtf(2 * 3.1415926f * sigma * sigma); float coef2 = -1.0 / (2 * sigma * sigma); float c_sum = 0; float f_sum = 0; int val0 = src[x + y * step]; for (int dy = 0; dy < 3; ++dy) { for (int dx = 0; dx < 3; ++dx) { int val = src[(x + dx) + (y + dy) * step]; int diff = val - val0; float w = filter3[dy][dx] * coef * expf(diff * diff * coef2); f_sum += w; c_sum += w * val; } } dst[x + y * step] = (int)(c_sum / f_sum + 0.5f); }} 必ず”f”を付ける。付けないと doubleの演 算になって、Tesla以外では、 かなり遅くなるので注意 重みの計算 21
22. Copyright© Fixstars Group バイラテラルフィルタ高速化 カーネルプロファイリング • かなり演算ネック 22
23. Copyright© Fixstars Group __global__ void BilateralKernelSimple( const uint8_t *src, uint8_t *dst, int width, int height, int step, float sigma) { int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; if (x < width && y < height) { float coef = 1.0 / sqrtf(2 * 3.1415926f * sigma * sigma); float coef2 = -1.0 / (2 * sigma * sigma); float c_sum = 0; float f_sum = 0; int val0 = src[x + y * step]; for (int dy = 0; dy < 3; ++dy) { for (int dx = 0; dx < 3; ++dx) { int val = src[(x + dx) + (y + dy) * step]; int diff = val - val0; float w = filter3[dy][dx] * coef * expf(diff * diff * coef2); f_sum += w; c_sum += w * val; } } dst[x + y * step] = (int)(c_sum / f_sum + 0.5f); }} 割り算、sqrt 割り算 exp 割り算 重い演算が多い 23
24. Copyright© Fixstars Group __global__ void BilateralKernelFast( const uint8_t *src, uint8_t *dst, int width, int height, int step, float sigma) { int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; if (x < width && y < height) { float coef = __frsqrt_rn(2 * 3.1415926f * sigma * sigma); float coef2 = __frcp_rn(-2 * sigma * sigma); float c_sum = 0; float f_sum = 0; int val0 = src[x + y * step]; for (int dy = 0; dy < 3; ++dy) { for (int dx = 0; dx < 3; ++dx) { int val = src[(x + dx) + (y + dy) * step]; int diff = val - val0; float w = filter3[dy][dx] * coef * __expf(diff * diff * coef2); f_sum += w; c_sum += w * val; } } dst[x + y * step] = (int)(__fdividef(c_sum, f_sum) + 0.5f); }} 1/sqrtf(x) のintrinsic 1/x のintrinsic expf(x) の高速版 x/y の高速版 24
25. Copyright© Fixstars Group バイラテラルフィルタ高速化 高速版の誤差 https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#intrinsic-functions ulp = Unit in the last place 25
26. Copyright© Fixstars Group バイラテラルフィルタ高速化 
 バイラテラルフィルタ 計算時間 (ms)
 最適化前
 2.48
 最適化後
 1.63
 最適化前
 最適化後
 26
27. Copyright© Fixstars Group 転置
28. Copyright© Fixstars Group 転置 28
29. Copyright© Fixstars Group 転置 単純に書いてみる __global__ void TransposeKernelSimple( const uint8_t *src, uint8_t *dst, int width, int height) { int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; if (x < width && y < height) dst[y + x * height] = src[x + y * width]; } xとyを逆にするだけ 29
30. Copyright© Fixstars Group 転置 カーネルプロファイル Loadに比べてStoreの L2⇔DRAM Transactionが かなり多い メモリがネック 30
31. Copyright© Fixstars Group 転置 単純に書いてみる __global__ void TransposeKernelSimple( const uint8_t *src, uint8_t *dst, int width, int height) { int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; if (x < width && y < height) dst[y + x * height] = src[x + y * width]; } 書き込みが全くコアレスアクセス になっていない 31
32. Copyright© Fixstars Group 転置 単純に書いてみる メモリ連続方向 メモリ連続方向 Warpのアクセスするメモリ コアレスアクセスになっている 完全に飛び飛びのアクセス Load Store NG！ 32
33. Copyright© Fixstars Group 転置 Shared Memoryを使う Warpのアクセスするメモリ Shared Memory Shared Memory 同じメモリ 書き込みもコアレスアクセス でできるようにする メモリ連続方向 メモリ連続方向 33
34. Copyright© Fixstars Group __global__ void TransposeKernelShared( const uint8_t *src, uint8_t *dst, int width, int height) { int tx = threadIdx.x; int ty = threadIdx.y; int xbase = blockIdx.x * blockDim.x; int ybase = blockIdx.y * blockDim.y; __shared__ uint8_t sbuf[16][16]; { int x = xbase + tx; int y = ybase + ty; if (x < width && y < height) sbuf[ty][tx] = src[x + y * width]; } __syncthreads(); { int x = xbase + ty; int y = ybase + tx; if (x < width && y < height) dst[y + x * height] = sbuf[tx][ty]; }} 転置 Shared Memoryを使う 書き込みもコアレスアクセスで できるようにする 一旦Shared Memoryに格納 34
35. Copyright© Fixstars Group 転置 Shared Memoryを使う LoadとStoreで L2⇔DRAM Transactionが 同じになった メモリネックでなくなった 35
36. Copyright© Fixstars Group 転置 Shared Memoryを使う 計測環境 CPU: Core i7-8700 3.2GHz (6コア 12スレッド) GPU: GeForce RTX 2060 計測条件 6720x4480の画像(グレースケール)を処理 計算時間のみで、データ転送やメモリ確保などの 時間を含めず 36
37. Copyright© Fixstars Group 転置 Shared Memoryを使う • ただし、Shared Memory のバンクコンフリクトが発生している 37
38. Copyright© Fixstars Group __global__ void TransposeKernelFast( const uint8_t *src, uint8_t *dst, int width, int height){ int tx = threadIdx.x; int ty = threadIdx.y; int xbase = blockIdx.x * blockDim.x; int ybase = blockIdx.y * blockDim.y; __shared__ uint8_t sbuf[16][16+4]; { int x = xbase + tx; int y = ybase + ty; if (x < width && y < height) sbuf[ty][tx] = src[x + y * width]; } __syncthreads(); { int x = xbase + ty; int y = ybase + tx; if (x < width && y < height) dst[y + x * height] = sbuf[tx][ty]; }} 転置 バンクコンフリクト回避 パディングを追加 Shared Memoryのバンクは 4バイトインターリーブされているので、 4バイトパディングを追加する 38
39. Copyright© Fixstars Group 転置 バンクコンフリクト回避 • バンクコンフリクトがゼロになった 39
40. Copyright© Fixstars Group 転置 バンクコンフリクト回避 計測環境 CPU: Core i7-8700 3.2GHz (6コア 12スレッド) GPU: GeForce RTX 2060 計測条件 6720x4480の画像(グレースケール)を処理 計算時間のみで、データ転送やメモリ確保などの 時間を含めず 40
41. Copyright© Fixstars Group 転置 1スレッドあたり処理量を増やす カーネルのアセンブラ • １スレッドあたりの処理量が少なすぎる ◦ カーネルが32命令しかない • １スレッドが処理する要素数を増やす
42. Copyright© Fixstars Group __global__ void TransposeKernelFast2( const uint8_t *src, uint8_t *dst, int width, int height){ int tx = threadIdx.x; int ty = threadIdx.y; int xbase = blockIdx.x * 32; int ybase = blockIdx.y * 32; __shared__ uint8_t sbuf[32][32+4]; { int x = xbase + tx; if (x < width) { int yend = min(ybase + 32, height); for (int tyy = ty, y = ybase + ty; y < yend; tyy += 8, y += 8) { sbuf[tyy][tx] = src[x + y * width]; }}} __syncthreads(); { int y = ybase + tx; if (y < height) { int xend = min(xbase + 32, width); for (int tyy = ty, x = xbase + ty; x < xend; tyy += 8, x += 8) { dst[y + x * height] = sbuf[tx][tyy]; }}}} 転置 1スレッドあたり処理量を増やす １スレッドが4要素処理するように修正 42
43. Copyright© Fixstars Group 転置 1スレッドあたり処理量を増やす 計測環境 CPU: Core i7-8700 3.2GHz (6コア 12スレッド) GPU: GeForce RTX 2060 計測条件 6720x4480の画像(グレースケール)を処理 計算時間のみで、データ転送やメモリ確保などの 時間を含めず 43
44. Copyright© Fixstars Group リダクション
45. Copyright© Fixstars Group リダクション • ２次元データのリダクション • ２つの軸のリダクションをそれぞれ実装してみる メモリ連続方向 x y Y軸リダクション X軸リダクション メモリ連続方向 45
46. Copyright© Fixstars Group リダクション Y軸リダクション • 1スレッド１列担当 • コアレスアクセスになっていることに注意 __global__ void ReduceHKernelSimple( const uint8_t *src, float *dst, int width, int height) { int x = blockIdx.x * blockDim.x + threadIdx.x; if (x < width) { float sum = 0; for (int y = 0; y < height; ++y) { sum += src[x + y * width]; } dst[x] = sum; } } Y軸リダクション 46
47. Copyright© Fixstars Group リダクション Y軸リダクション • 6720x4480だとブロック数が7個しかない ◦ 1スレッド１列担当なので ceil(6720/1024)=7 47
48. Copyright© Fixstars Group リダクション Y軸リダクション • 列を分割して並列数を増やす ◦ 1列1スレッド→ceil(行数/128)スレッド __global__ void ReduceHKernelFast( const uint8_t *src, float *dst, int width, int height) { int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * 128; if (x < width) { float sum = 0; for (int yend = min(y + 128, height); y < yend; ++y) { sum += src[x + y * width]; } atomicAdd(&dst[x], sum); } } 128行ごとに分割して処理する dstはこのカーネルを呼び出す前に ゼロ初期化しておく 48
49. Copyright© Fixstars Group リダクション Y軸リダクション 計測環境 CPU: Core i7-8700 3.2GHz (6コア 12スレッド) GPU: GeForce RTX 2060 計測条件 6720x4480の画像(グレースケール)を処理 計算時間のみで、データ転送やメモリ確保などの 時間を含めず 49
50. Copyright© Fixstars Group リダクション X軸リダクション • Y軸リダクションと同じように実装 X軸リダクション __global__ void ReduceWKernelSimple( const uint8_t *src, float *dst, int width, int height) { int y = blockIdx.x * blockDim.x + threadIdx.x; int x = blockIdx.y * 128; if (y < height) { float sum = 0; for (int xend = min(x + 128, width); x < xend; ++x) { sum += src[x + y * width]; } atomicAdd(&dst[y], sum); } } 50
51. Copyright© Fixstars Group リダクション X軸リダクション • Y軸リダクションと同じように実装 • 2.39ms …遅い X軸リダクション 
 X軸リダクション時間 (ms)
 Y軸と同じ方法
 2.39
 
 Y軸リダクション時間 (ms)
 1列1スレッド
 0.28
 列も分割
 0.115
 参考 51
52. Copyright© Fixstars Group リダクション X軸リダクション __global__ void ReduceWKernelSimple( const uint8_t *src, float *dst, int width, int height) { int y = blockIdx.x * blockDim.x + threadIdx.x; int x = blockIdx.y * 128; if (y < height) { float sum = 0; for (int xend = min(x + 128, width); x < xend; ++x) { sum += src[x + y * width]; } atomicAdd(&dst[y], sum); } } このアクセスが全く コアレスアクセスでない 52
53. Copyright© Fixstars Group リダクション X軸リダクション X軸リダクション このアクセスはダメ Warpのアクセスするメモリ こうアクセスしたい！ 水平方向の リダクション が必要 53
54. Copyright© Fixstars Group リダクション 水平方向のリダクションが必要 パラレルリダクション 54
55. Copyright© Fixstars Group リダクション パラレルリダクション • 1行を1ブロックが担当 __global__ void ReduceWKernelFast( const uint8_t *src, float *dst, int width, int height) { int tid = threadIdx.x; int y = blockIdx.y; __shared__ float sbuf[512]; float sum = 0; for (int x = tid; x < width; x += 512) { sum += src[x + y * width]; } sbuf[tid] = sum; __syncthreads(); sum = ReduceFunc(tid, sbuf); if (tid == 0) dst[y] = sum; } 512要素までのリダクションは普通にス レッドごとに計算 1ブロック512スレッドで コードを書いた場合 Shared Memoryに書いて パラレルリダクションを呼び出す 55
56. Copyright© Fixstars Group リダクション パラレルリダクション __device__ float ReduceFunc(int tid, float* buf) { if (tid < 256) { buf[tid] += buf[tid + 256]; } __syncthreads(); if (tid < 128) { buf[tid] += buf[tid + 128]; } __syncthreads(); if (tid < 64) { buf[tid] += buf[tid + 64]; } __syncthreads(); float sum; if (tid < 32) { sum = buf[tid] + buf[tid + 32]; sum += __shfl_down_sync(0xffffffff, sum, 16); sum += __shfl_down_sync(0xffffffff, sum, 8); sum += __shfl_down_sync(0xffffffff, sum, 4); sum += __shfl_down_sync(0xffffffff, sum, 2); sum += __shfl_down_sync(0xffffffff, sum, 1); } return sum; } 32スレッドまでは __syncthreads()を使って計算 32スレッドになったら、Warp Shuffleで計算 56
57. Copyright© Fixstars Group リダクション パラレルリダクション 57
58. Copyright© Fixstars Group 本セミナーのまとめ • バイラテラルフィルタ高速化 ◦ CUDA組み込み関数を使って演算を軽量化 • 転置 ◦ Shared Memoryを使ったメモリアクセス最適化 ◦ バンクコンフリクト回避 • リダクション ◦ X軸、Y軸方向のリダクション ◦ メモリのアクセス方向を意識した計算 ◦ 水平方向のリダクションを高速に行うパラレルリダクション 58
59. Copyright © Fixstars Group Thank you! お問い合わせ窓口 : hr-seminar@fixstars.com