CUDA入門 #2

この記事は書籍「Hands-On GPU Programming with Python and CUDA」学習時に生じた疑問点などをChatGPT 4oに質問した内容をもとにChatGPT 4oにブログ記事としてまとめてもらったものです

1. PyCUDAのgpuarrayクラスの概要

CUDAを使用する際、ホスト（CPU）からGPUにデータを転送し、GPU上で処理を行う必要があります。PyCUDAは、この操作を簡単に行うためのPythonラッパーを提供しており、その中でもgpuarrayクラスは非常に重要な役割を担っています。gpuarrayクラスは、NumPy配列と似たインターフェースを持ち、GPU上で配列操作を簡単に行えるように設計されています。

基本的な操作

gpuarrayは、NumPy配列をGPUメモリに転送するための便利なツールです。例えば、以下のようにNumPy配列をgpuarrayに変換し、GPU上で操作できます。

import pycuda.autoinit
import pycuda.gpuarray as gpuarray
import numpy as np

a = np.random.randn(5).astype(np.float32)  # NumPy配列の作成
a_gpu = gpuarray.to_gpu(a)                 # NumPy配列をGPUに転送

GPU上にデータが転送された後、通常の演算（加算や乗算など）はGPU上で直接行われます。例えば、次のコードではGPU上の2つの配列を加算しています。

b_gpu = gpuarray.to_gpu(np.random.randn(5).astype(np.float32))
c_gpu = a_gpu + b_gpu  # GPU上での配列演算

gpuarrayを使えば、メモリ管理を手動で行うことなく、非常にシンプルな形でGPU上でのデータ処理を行うことができます。

結果の取得

GPU上で計算された結果をホストに戻す場合、get()メソッドを使用します。これにより、GPU上のデータをホスト側のNumPy配列として取得できます。

c = c_gpu.get()  # GPUからCPUへデータを転送

高度な数値演算

gpuarrayは基本的な四則演算だけでなく、dot（内積）、exp（指数関数）、sin（正弦関数）などの数学関数もサポートしています。これにより、複雑な数値演算を簡単に実装することが可能です。

このように、gpuarrayクラスは、GPU上での効率的な数値演算を簡潔に実装するための便利なツールです。続いて、gpuarrayとカスタムカーネルの連携について紹介します。

データ型

NumPy dtype	C/C++の型	説明
np.int8	char	8ビット符号付き整数
np.uint8	unsigned char	8ビット符号なし整数
np.int16	short	16ビット符号付き整数
np.uint16	unsigned short	16ビット符号なし整数
np.int32	int	32ビット符号付き整数
np.uint32	unsigned int	32ビット符号なし整数
np.int64	long long	64ビット符号付き整数
np.uint64	unsigned long long	64ビット符号なし整数
np.float32	float	単精度浮動小数点（32ビット）
np.float64	double	倍精度浮動小数点（64ビット）
np.bool_	bool	真偽値

gpuarrayはCUDAの標準的なメモリ管理と同じく、ホストメモリとは分離されていますが、cudaMallocManagedのような統一メモリを使う場合には直接的な管理が不要です。どちらを使うかはユースケースによります。

2. カスタムカーネルとの連携

gpuarrayは、シンプルな数値演算に便利なツールですが、より複雑な処理を行いたい場合にはカスタムカーネルとの連携が不可欠です。カスタムカーネルは、CUDAコードで記述され、GPU上で並列実行される関数です。PyCUDAでは、これを簡単に統合して使用することができます。

カスタムカーネルの定義

PyCUDAのSourceModuleを使って、カスタムCUDAカーネルを定義し、Pythonから呼び出すことができます。次に、簡単な例として、配列の各要素に定数を掛けるカーネルを定義して実行する方法を紹介します。

import pycuda.autoinit
import pycuda.driver as cuda
import pycuda.gpuarray as gpuarray
import numpy as np
from pycuda.compiler import SourceModule

# CUDAカーネルの定義
mod = SourceModule("""
__global__ void multiply_by_scalar(float *a, float scalar)
{
    int idx = threadIdx.x;
    a[idx] = a[idx] * scalar;
}
""")

# カーネル関数を取得
multiply_by_scalar = mod.get_function("multiply_by_scalar")

# GPU上で計算を行うためのデータを準備
a = np.array([1, 2, 3, 4, 5], dtype=np.float32)
a_gpu = gpuarray.to_gpu(a)

# カーネルを呼び出す
scalar_value = np.float32(10.0)
multiply_by_scalar(a_gpu, scalar_value, block=(a_gpu.size, 1, 1))

# 計算結果を取得
result = a_gpu.get()
print("計算結果:", result)

この例では、配列aの各要素に対して定数10を掛けています。SourceModuleを使うことで、C言語風のCUDAカーネルをPython内で定義し、gpuarrayと連携してデータのやり取りが可能です。

カーネル実行のポイント

CUDAカーネルの定義: __global__修飾子を使ってカーネルを定義します。これは、GPU上で実行される関数であることを示します。
スレッドとブロックの指定: カーネル呼び出し時にblock=(スレッド数, 1, 1)のようにブロック内のスレッド数を指定します。この設定により、どのように並列処理が行われるかが決まります。
データのやり取り: gpuarray.to_gpu()でホストからGPUにデータを転送し、計算結果をgpuarray.get()でホスト側に戻すことができます。

カスタムカーネルを使うことで、gpuarrayを超えた高度な並列処理が可能になります。次は、要素ごとの演算を簡単に記述できるElementwiseKernelについて紹介します。

3. ElementwiseKernelの使い方と活用方法

ElementwiseKernelは、PyCUDAが提供する非常に便利なツールで、カーネルコードをシンプルに記述し、GPU上での要素ごとの演算を効率的に実行できる機能です。カーネルコードを自分で詳細に記述する必要がないため、複雑なカーネルを書かずとも、各要素に対して同じ演算を行うような処理を簡潔に実装できます。

ElementwiseKernelの基本構造

まず、ElementwiseKernelの基本構造を見てみましょう。次のように、Python内で要素ごとの処理を記述し、並列に実行できます。

from pycuda.elementwise import ElementwiseKernel
import pycuda.gpuarray as gpuarray
import numpy as np

# 2つの配列の各要素を加算するカーネル
add_kernel = ElementwiseKernel(
    "float *a, float *b, float *c",  # 3つの引数（2つの入力と1つの出力）
    "c[i] = a[i] + b[i];",           # 要素ごとの加算を定義
    "add_arrays"                     # カーネル名
)

# NumPy配列を準備し、GPUに転送
a = np.array([1, 2, 3, 4, 5], dtype=np.float32)
b = np.array([10, 20, 30, 40, 50], dtype=np.float32)
a_gpu = gpuarray.to_gpu(a)
b_gpu = gpuarray.to_gpu(b)
c_gpu = gpuarray.empty_like(a_gpu)  # 結果を格納する配列

# カーネルを実行
add_kernel(a_gpu, b_gpu, c_gpu)

# 結果をホストに転送して表示
result = c_gpu.get()
print("加算結果:", result)

このコードでは、a[i]とb[i]を加算してc[i]に格納するというシンプルな要素ごとの処理を行っています。このElementwiseKernelを使うことで、冗長なカーネルコードを書くことなく、並列処理を簡単に実現できます。

ElementwiseKernelのメリット

簡潔な記述: カーネルコードをシンプルに記述でき、わざわざSourceModuleを使ってカーネル関数を詳細に定義する必要がありません。

• 自動スレッド管理: ElementwiseKernelは自動的にスレッドとブロックを最適化してくれるため、手動で細かくスレッド管理をする必要がなくなります。
• 柔軟性: 基本的な加算や乗算だけでなく、複雑な条件付き演算や数学関数を使った演算も実装可能です。
• 応用例: 条件付き処理

次に、もう少し複雑な処理として、条件に基づいて演算を行う例を見てみましょう。以下では、各要素の値が5より大きい場合に10を加えるという処理を行います。

complex_kernel = ElementwiseKernel(
    "float *a, float *result",  # 引数を定義
    """
    float temp = a[i];
    if (temp > 5.0) {
        result[i] = temp + 10;
    } else {
        result[i] = temp;
    }
    """,
    "complex_kernel"
)

# NumPy配列を準備
a = np.array([1, 4, 6, 9, 2], dtype=np.float32)
a_gpu = gpuarray.to_gpu(a)
result_gpu = gpuarray.empty_like(a_gpu)

# カーネルを実行
complex_kernel(a_gpu, result_gpu)

# 結果を取得して表示
result = result_gpu.get()
print("条件付き演算結果:", result)

このように、ElementwiseKernelを使えば、条件分岐を含む複雑な処理も簡単に並列化できます。

ElementwiseKernelは、要素ごとの処理を簡単に記述し、GPUの並列処理の力を効率よく活用するための強力なツールです。次は、配列の累積和などの処理に便利なInclusiveScanKernelについて紹介します。

4. InclusiveScanKernelの紹介と例

InclusiveScanKernelは、配列内の要素を累積していく操作を効率的に実行するためのカーネルジェネレータです。スキャン操作は、各要素に対してその位置までのすべての要素を含めた累積結果を計算する操作で、CUDAの並列計算で頻繁に使用されます。特に累積和や累積積の計算に便利です。

Inclusive Scanとは？

Inclusive Scanは、現在の要素までの累積和を計算する操作です。例えば、次の配列に対して：

[1, 2, 3, 4]

Inclusive Scanの結果は次のようになります：

[1, 3, 6, 10]

これはそれぞれ、1、1+2、1+2+3、1+2+3+4 という累積結果です。

InclusiveScanKernelの基本構造
PyCUDAでは、この累積操作を簡単に実装するためにInclusiveScanKernelが提供されています。以下が、その基本的な使用方法です。

from pycuda.scan import InclusiveScanKernel
import pycuda.gpuarray as gpuarray
import numpy as np

# InclusiveScanKernelを定義（加算による累積和）
scan_kernel = InclusiveScanKernel(np.float32, "a + b")

# NumPy配列を用意し、GPUに転送
a = np.array([1, 2, 3, 4, 5], dtype=np.float32)
a_gpu = gpuarray.to_gpu(a)

# カーネルを実行
scan_kernel(a_gpu)

# 結果をホストに転送して表示
result = a_gpu.get()
print("累積和:", result)

出力結果:

累積和: [ 1.  3.  6. 10. 15.]

この例では、配列の各要素に対して累積和が計算されています。InclusiveScanKernelを使えば、数行のコードで高速なスキャン操作を実装することが可能です。

他の演算の例

InclusiveScanKernelは、累積和だけでなく、累積積や最大値、最小値など様々な演算に応用できます。次に、累積積を計算する例を示します。

# 累積積を計算するカーネル
scan_kernel = InclusiveScanKernel(np.float32, "a * b")

# NumPy配列を用意し、GPUに転送
a = np.array([1, 2, 3, 4, 5], dtype=np.float32)
a_gpu = gpuarray.to_gpu(a)

# カーネルを実行
scan_kernel(a_gpu)

# 結果をホストに転送して表示
result = a_gpu.get()
print("累積積:", result)

出力結果:

累積積: [ 1.  2.  6. 24. 120.]

この例では、配列の要素ごとの累積積が計算されます。InclusiveScanKernelは、非常に簡単にカスタムスキャン操作を行えるため、様々な用途に応用できます。

まとめ

InclusiveScanKernelは、累積演算を効率的に並列化するための強力なツールです。累積和や累積積といった操作はもちろんのこと、他の演算も柔軟に実装できるため、特に大規模データの処理や並列化された計算において非常に役立ちます。

次は、他の関連カーネル（ExclusiveScanKernelやReductionKernelなど）について紹介します。

5. CUDAの他の関連カーネル

CUDAには、InclusiveScanKernel以外にも、様々な演算を効率的に実行するための便利なカーネルが用意されています。これらのカーネルを活用することで、複雑な演算を高速かつ簡単に並列化することができます。ここでは、いくつかの代表的なカーネルについて紹介します。

1. ExclusiveScanKernel

ExclusiveScanKernelは、InclusiveScanKernelと同様のスキャン操作を行いますが、現在の要素を除いた累積結果を計算するカーネルです。これを「除外スキャン（exclusive scan）」と呼びます。たとえば、以下のような配列に対して：

[1, 2, 3, 4]

Exclusive Scanの結果は：

[0, 1, 3, 6]

例: ExclusiveScanKernelによる累積和の計算

from pycuda.scan import ExclusiveScanKernel
import pycuda.gpuarray as gpuarray
import numpy as np

# ExclusiveScanKernelを定義（加算による累積和）
exclusive_scan_kernel = ExclusiveScanKernel(np.float32, "a + b")

# NumPy配列を準備し、GPUに転送
a = np.array([1, 2, 3, 4, 5], dtype=np.float32)
a_gpu = gpuarray.to_gpu(a)

# カーネルを実行
exclusive_scan_kernel(a_gpu)

# 結果をホストに転送して表示
result = a_gpu.get()
print("Exclusive累積和:", result)

出力結果:

Exclusive累積和: [ 0.  1.  3.  6. 10.]

ExclusiveScanKernelは、特定の処理に対して「前の要素までの結果が必要」な場合に便利です。

2. ReductionKernel

ReductionKernelは、配列全体を1つの結果に縮約（reduce）するためのカーネルです。これにより、配列の合計や最大値、最小値などを並列に計算できます。

例: 配列の合計を計算するReductionKernel

from pycuda.reduction import ReductionKernel
import pycuda.gpuarray as gpuarray
import numpy as np

# ReductionKernelを定義（加算による縮約）
reduction_kernel = ReductionKernel(
    np.float32, neutral="0",
    reduce_expr="a+b", map_expr="x[i]",
    arguments="float *x"
)

# NumPy配列を準備し、GPUに転送
a = np.random.randn(10).astype(np.float32)
a_gpu = gpuarray.to_gpu(a)

# カーネルを実行して配列の合計を計算
result = reduction_kernel(a_gpu).get()

# 結果を表示
print("配列の合計:", result)

出力結果:

配列の合計: [配列の全要素の合計]

ReductionKernelは、加算以外にも最大値、最小値、乗算などの演算を行うことが可能です。

3. SegmentedScanKernel

SegmentedScanKernelは、セグメント（区切り）ごとにスキャン操作を行うカーネルです。通常のスキャン操作をセグメントごとに分けて実行する際に使用します。

例: セグメントごとの累積和

from pycuda.scan import SegmentedScanKernel
import pycuda.gpuarray as gpuarray
import numpy as np

# NumPy配列とセグメントを準備
a = np.array([1, 2, 3, 4, 5, 1, 2, 3, 4, 5], dtype=np.float32)
segment_flags = np.array([1, 0, 0, 0, 0, 1, 0, 0, 0, 0], dtype=np.int32)

# SegmentedScanKernelを定義（加算による累積和）
segmented_scan_kernel = SegmentedScanKernel(np.float32, np.int32, "a + b")

# GPUにデータを転送
a_gpu = gpuarray.to_gpu(a)
segment_flags_gpu = gpuarray.to_gpu(segment_flags)

# カーネルを実行
segmented_scan_kernel(a_gpu, segment_flags_gpu)

# 結果をホストに転送して表示
result = a_gpu.get()
print("セグメントごとの累積和:", result)

出力結果:

セグメントごとの累積和: [ 1.  3.  6. 10. 15.  1.  3.  6. 10. 15.]

SegmentedScanKernelを使うと、特定のセグメントごとに演算を行いたい場合に非常に便利です。

まとめ

InclusiveScanKernel以外にも、CUDAには様々なカーネルが用意されており、異なる演算に対して効率的に並列処理を実行できます。ExclusiveScanKernelで除外スキャンを行ったり、ReductionKernelで配列全体を1つの結果にまとめたりと、さまざまな操作を簡単に実装できます。次は、画像処理向けのカーネルについて詳しく見ていきます。

6. 画像処理向けのCUDAカーネル

CUDAは、並列処理が得意なため、画像処理に非常に適しています。特に、画像に対してフィルタ処理やピクセル単位の操作を行う場合、CUDAカーネルを利用することで大幅に処理速度を向上させることができます。ここでは、画像処理に向いているカーネルや、その活用方法について紹介します。

1. ElementwiseKernelを使った画像処理

画像処理の多くは、各ピクセルに対して同じ操作を行う場合が多いため、ElementwiseKernelを使用することでシンプルに並列処理が可能です。たとえば、画像の反転処理や輝度の調整などの単純な操作は、ElementwiseKernelで簡単に実装できます。

例: 画像の反転処理

import pycuda.autoinit
import pycuda.gpuarray as gpuarray
from pycuda.elementwise import ElementwiseKernel
import numpy as np

# 256x256ピクセルのランダム画像を生成
image = np.random.randint(0, 256, (256, 256)).astype(np.float32)

# GPUに転送
image_gpu = gpuarray.to_gpu(image)

# ElementwiseKernelを使用して画像の反転処理を実装
invert_kernel = ElementwiseKernel(
    "float *image",              # 1つの入力（画像）
    "image[i] = 255 - image[i];", # 各ピクセルの値を反転
    "invert_image"               # カーネル名
)

# カーネルを実行
invert_kernel(image_gpu)

# 結果をホストに転送して表示
inverted_image = image_gpu.get()
print("反転処理後の画像:", inverted_image)

このコードでは、256×256ピクセルの画像の各ピクセルの値を反転しています。ElementwiseKernelは要素ごとの操作に非常に適しており、ピクセル単位のシンプルな画像処理に活用できます。

2. カスタムカーネルを使った畳み込みフィルタ

画像処理では、畳み込みフィルタを使ったスムージングやエッジ検出がよく行われます。畳み込みフィルタでは、各ピクセルに対してその周囲のピクセルを使って計算を行います。これはCUDAカーネルを使って並列化することで、効率的に実行できます。

例: 3x3の畳み込みフィルタ
次に、CUDAカーネルを使って3x3の畳み込みフィルタを実装します。このフィルタは、各ピクセルとその周囲のピクセルを使ってスムージング（ぼかし）処理を行います。

import pycuda.autoinit
import pycuda.driver as cuda
import pycuda.gpuarray as gpuarray
import numpy as np
from pycuda.compiler import SourceModule

# 256x256ピクセルのランダム画像を生成
image = np.random.randint(0, 256, (256, 256)).astype(np.float32)

# 畳み込みカーネル (3x3のスムージングフィルタ)
kernel = np.array([[1/9, 1/9, 1/9], 
                   [1/9, 1/9, 1/9], 
                   [1/9, 1/9, 1/9]], dtype=np.float32)

# GPUに画像とカーネルを転送
image_gpu = gpuarray.to_gpu(image)
kernel_gpu = gpuarray.to_gpu(kernel)

# CUDAカーネル定義
mod = SourceModule("""
__global__ void convolve(float *image, float *output, float *kernel, int width, int height)
{
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    
    if (x >= width || y >= height)
        return;

    float sum = 0.0;
    int kernel_size = 3;
    int half_size = kernel_size / 2;

    // 畳み込み演算
    for (int i = -half_size; i <= half_size; i++) {
        for (int j = -half_size; j <= half_size; j++) {
            int ix = min(max(x + i, 0), width - 1);
            int iy = min(max(y + j, 0), height - 1);
            sum += image[iy * width + ix] * kernel[(i + half_size) * kernel_size + (j + half_size)];
        }
    }
    output[y * width + x] = sum;
}
""")

# カーネル関数を取得
convolve = mod.get_function("convolve")

# 結果を格納する出力配列をGPU上に用意
output_gpu = gpuarray.empty_like(image_gpu)

# カーネルを実行
block_size = (16, 16, 1)
grid_size = (16, 16, 1)
convolve(image_gpu, output_gpu, kernel_gpu, np.int32(image.shape[1]), np.int32(image.shape[0]), block=block_size, grid=grid_size)

# 結果をホストに転送
output_image = output_gpu.get()
print("畳み込み後の画像:", output_image)

この例では、3×3のスムージングフィルタを画像に適用しています。各スレッドが1つのピクセルを担当し、その周辺のピクセルを使って計算を行います。

3. 共有メモリの活用

CUDAには**共有メモリ（shared memory）**という機能があり、これを活用することで畳み込みフィルタのようなローカルなデータアクセスを高速化することができます。スレッドブロック内のスレッド同士でデータを共有することで、メモリアクセスの回数を減らし、より効率的な処理が可能になります。

例えば、畳み込みフィルタを実行する際に、共有メモリを使用して各スレッドが担当する領域のデータを一度共有メモリに読み込むことで、重複するメモリアクセスを最小化できます。

まとめ

画像処理において、CUDAの並列処理は非常に強力です。ピクセル単位の操作を効率化するElementwiseKernelや、畳み込みフィルタのような複雑な操作もカスタムカーネルで実装することが可能です。また、共有メモリを使うことで、さらに高度な最適化を行うことができ、画像処理のパフォーマンスを大幅に向上させることができます。

次は、カーネルの定義と実行に便利なSourceModuleの使い方と応用について紹介します。

7. SourceModuleの使い方と応用

CUDAカーネルを使用する際に、PyCUDAではSourceModuleが非常に重要な役割を果たします。SourceModuleを使用することで、C言語風のCUDAカーネルコードをPython内で直接定義し、GPU上で実行することができます。これにより、カーネルコードのコンパイルと実行をスムーズに行うことができ、複雑な計算処理を効率よく実装できます。

SourceModuleの基本構造
SourceModuleを使うことで、以下のようにCUDAカーネルをPythonコードの中に記述し、実行できます。

from pycuda.compiler import SourceModule
import pycuda.autoinit
import pycuda.driver as cuda
import numpy as np

# CUDAカーネルの定義
mod = SourceModule("""
__global__ void add_constant(float *a, float constant)
{
    int idx = threadIdx.x;
    a[idx] += constant;
}
""")

# カーネル関数を取得
add_constant = mod.get_function("add_constant")

# NumPy配列を作成
a = np.random.randn(10).astype(np.float32)

# GPUメモリに転送
a_gpu = cuda.mem_alloc(a.nbytes)
cuda.memcpy_htod(a_gpu, a)

# カーネルを実行
constant_value = np.float32(5.0)
add_constant(a_gpu, constant_value, block=(10, 1, 1))

# 結果を取得して表示
cuda.memcpy_dtoh(a, a_gpu)
print("結果:", a)

この例では、各スレッドが配列aの要素に対して定数5.0を加算するカーネルを実行しています。

SourceModuleのポイント
カーネルの定義： カーネルコードをSourceModuleの引数としてC言語風に記述します。この場合、__global__修飾子を付けてカーネル関数を定義しています。この関数は、GPU上で実行され、各スレッドが並行して同じ処理を行います。

カーネル関数の取得： mod.get_function("カーネル名")を使って、定義したカーネルを取得します。この関数を使用して、Pythonコード内でカーネルを呼び出し、引数を渡して実行できます。

スレッドとブロックの指定： カーネル呼び出し時に、block=(スレッド数, 1, 1)といった形式でスレッド数やブロックサイズを指定します。これにより、カーネルの並列処理の仕方を指定できます。

複数の引数を持つカーネル
次に、複数の配列を操作するカーネルを実装する例を見てみましょう。この例では、2つの配列の各要素を加算し、結果を別の配列に格納します。

# CUDAカーネルの定義
mod = SourceModule("""
__global__ void add_arrays(float *a, float *b, float *c)
{
    int idx = threadIdx.x;
    c[idx] = a[idx] + b[idx];
}
""")

# カーネル関数を取得
add_arrays = mod.get_function("add_arrays")

# NumPy配列を準備
a = np.array([1, 2, 3, 4, 5], dtype=np.float32)
b = np.array([10, 20, 30, 40, 50], dtype=np.float32)

# GPUメモリに転送
a_gpu = cuda.mem_alloc(a.nbytes)
b_gpu = cuda.mem_alloc(b.nbytes)
c_gpu = cuda.mem_alloc(a.nbytes)

cuda.memcpy_htod(a_gpu, a)
cuda.memcpy_htod(b_gpu, b)

# カーネルを実行
add_arrays(a_gpu, b_gpu, c_gpu, block=(a.size, 1, 1))

# 結果を取得
c = np.empty_like(a)
cuda.memcpy_dtoh(c, c_gpu)

print("加算結果:", c)

このコードでは、2つの配列aとbの各要素を加算し、その結果を配列cに格納しています。

応用例: 行列乗算

次に、より複雑な例として行列の乗算を行うカーネルを実装してみます。行列乗算は、並列処理に適しており、CUDAを使うと効率的に計算が行えます。

# CUDAカーネルの定義
mod = SourceModule("""
__global__ void matrix_multiply(float *a, float *b, float *c, int width)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    
    float sum = 0.0;
    for (int k = 0; k < width; ++k)
    {
        sum += a[row * width + k] * b[k * width + col];
    }
    c[row * width + col] = sum;
}
""")

# カーネル関数を取得
matrix_multiply = mod.get_function("matrix_multiply")

# 4x4の行列を準備
a = np.random.randn(4, 4).astype(np.float32)
b = np.random.randn(4, 4).astype(np.float32)

# GPUメモリに転送
a_gpu = cuda.mem_alloc(a.nbytes)
b_gpu = cuda.mem_alloc(b.nbytes)
c_gpu = cuda.mem_alloc(a.nbytes)

cuda.memcpy_htod(a_gpu, a)
cuda.memcpy_htod(b_gpu, b)

# カーネルを実行
block_size = (2, 2, 1)
grid_size = (2, 2, 1)
matrix_multiply(a_gpu, b_gpu, c_gpu, np.int32(a.shape[1]), block=block_size, grid=grid_size)

# 結果を取得して表示
c = np.empty_like(a)
cuda.memcpy_dtoh(c, c_gpu)

print("行列乗算結果:", c)

この例では、2つの4×4行列の乗算を行っています。各スレッドが行列の特定の位置を担当し、部分的な計算を行うことで効率的に処理を進めています。

まとめ

SourceModuleを使えば、CUDAカーネルをPythonコード内で簡単に定義し、強力な並列計算を実行できます。シンプルな配列操作から複雑な行列演算まで、幅広い処理に応用できるため、効率的な並列処理を実現するための重要なツールです。

次は、CUDAの基本概念であるスレッド、ブロック、グリッドについて詳しく説明します。

8. CUDAのスレッド、ブロック、グリッドの説明

CUDAプログラミングの基本的な概念として、スレッド（Thread）、ブロック（Block）、**グリッド（Grid）**があります。これらは、GPU上で並列処理を行う際の処理単位です。それぞれがどのように機能し、どのように協調して並列計算を実行するかを理解することは、CUDAプログラムの設計において非常に重要です。ここでは、工場の作業に例えて分かりやすく説明します。

1. スレッド（Thread）＝ 作業員

スレッドは、CUDAにおける最も小さな処理単位で、1人の作業員に例えられます。スレッドは、特定のタスクを担当し、独立して計算を行います。例えば、画像処理では、スレッドが1つのピクセルに対する処理を担当します。

例：工場で100個の部品を組み立てる場合、100人の作業員がそれぞれ1個ずつ担当して組み立てを行います。同様に、CUDAでは100個のデータを100のスレッドが並列に処理します。

2. ブロック（Block）＝ 作業チーム

ブロックは、複数のスレッドをまとめた単位であり、作業チームに相当します。同じブロック内のスレッドは協力して作業を進めることができ、ブロック内ではデータの共有や同期が可能です。

例：10人の作業員が1つの大型製品の異なるパーツを組み立てるチームを考えてください。このチームは1つのブロックに相当し、各作業員（スレッド）はそのチーム内で協力して製品を組み立てます。
ブロック内のスレッドは、**共有メモリ（shared memory）**を使用してデータをやり取りすることができます。

3. グリッド（Grid）＝ 複数の作業チームからなる工場全体

グリッドは、複数のブロックをまとめた単位で、工場全体に相当します。複数の作業チームが同時に作業を行い、それぞれ異なる製品や異なるパーツを担当します。同様に、CUDAのグリッドは、全体の処理を並行して進めます。

例：工場には10チームがあり、それぞれのチームが10個の製品を組み立てるとします。全体として100個の製品が並行して組み立てられます。

グリッド全体で大量のデータを処理するために、複数のブロックが協力して作業を行います。
スレッド、ブロック、グリッドの関係
CUDAプログラムでは、各スレッドが個々のデータ（例えば、画像の1ピクセル）を処理し、同じブロック内のスレッドは連携して作業します。さらに、複数のブロックが協力して、全体の大規模データセットに対して並列処理を行います。

具体例：画像処理の工場に例える

256×256ピクセルの画像処理を考えます。

• スレッド：各スレッドは1ピクセルの処理を担当。
• ブロック：1つのブロックが16×16ピクセルの範囲を担当し、その範囲内の256個のピクセルを処理する。
• グリッド：16×16ブロックで構成され、全体で256×256ピクセルの画像を処理する。

このように、スレッド、ブロック、グリッドが階層的に構成され、効率的に並列処理を行います。

まとめ

CUDAのスレッド、ブロック、グリッドは、並列処理を効率化するための基本的な構成要素です。スレッドは個々の作業員としてデータを処理し、ブロックは作業チームとして協力しながら処理を進め、グリッドは工場全体として大規模なデータセットに対して並列処理を実行します。これらの概念をうまく活用することで、非常に効率的な並列計算を実現することができます。

次は、スレッド間での同期処理に関する重要な機能である__syncthreads()について紹介します。

9. スレッド間の同期：__syncthreads()の活用

CUDAプログラムにおける並列処理では、スレッドが独立して並行に動作しますが、時にはスレッド間で協調して作業を進める必要があります。このときに重要な役割を果たすのが__syncthreads()という関数です。__syncthreads()は、同じブロック内の全てのスレッドが同じ地点に到達するまで待機し、同期を取るための機能です。これにより、全てのスレッドが次のステップに進む前に、正しくデータが共有されたり、計算が完了したりすることが保証されます。

__syncthreads()の役割
__syncthreads()は、以下の2つの重要な役割を担っています：

メモリの同期：スレッドが共有メモリに書き込んだデータが全てのスレッドで整合性を保った状態で読み取れるようにします。
実行の同期：全てのスレッドが同じ箇所に到達するまで、次の処理を開始しないようにします。
これにより、データの不整合や競合状態を防ぎ、正しい結果が得られるようになります。

例：共有メモリを使った加算処理
次に、__syncthreads()を使って、スレッド間のデータを正しく同期しながら処理を進める例を紹介します。以下の例では、各スレッドがグローバルメモリから共有メモリにデータをコピーし、そのデータを使用して加算処理を行っています。

import pycuda.autoinit
import pycuda.driver as cuda
import numpy as np
from pycuda.compiler import SourceModule

# CUDAカーネルの定義
mod = SourceModule("""
__global__ void add_shared_memory(float *input, float *output)
{
    __shared__ float shared_data[10];  // 共有メモリの宣言
    int idx = threadIdx.x;

    // グローバルメモリから共有メモリにデータをコピー
    shared_data[idx] = input[idx];

    // 全スレッドが共有メモリへの書き込みを完了するのを待つ
    __syncthreads();

    // 共有メモリのデータを使って計算
    output[idx] = shared_data[idx] + 10;
}
""")

# カーネル関数を取得
add_shared_memory = mod.get_function("add_shared_memory")

# NumPy配列を準備
input_data = np.random.randn(10).astype(np.float32)

# GPUメモリに転送
input_gpu = cuda.mem_alloc(input_data.nbytes)
cuda.memcpy_htod(input_gpu, input_data)

# 結果を格納するための配列をGPUに用意
output_gpu = cuda.mem_alloc(input_data.nbytes)

# カーネルを実行
add_shared_memory(input_gpu, output_gpu, block=(10, 1, 1))

# 結果をホストに転送して表示
output_data = np.empty_like(input_data)
cuda.memcpy_dtoh(output_data, output_gpu)

print("入力データ:", input_data)
print("結果データ:", output_data)

例の解説

• 共有メモリの宣言： カーネル内で__shared__キーワードを使い、共有メモリ領域を宣言しています。このメモリはブロック内の全てのスレッドで共有され、効率的にデータをやり取りすることができます。
• 同期の実行： 各スレッドが共有メモリにデータを書き込んだ後、__syncthreads()を使って全スレッドが書き込みを終えるのを待機します。これにより、すべてのスレッドが確実に共有メモリにデータを保存してから、次のステップでそのデータを使って計算を行います。
• 同期後の処理： __syncthreads()で同期が完了した後、各スレッドは共有メモリのデータを使って計算を続行します。この例では、各要素に10を加えています。

__syncthreads()を使う場面

• 共有メモリを使う場合： 共有メモリを使用して、複数のスレッドが同じデータにアクセス・操作する際には、__syncthreads()で全スレッドを同期させることが必要です。
• 競合状態を防ぐため： スレッドが同じメモリ領域に対して異なるタイミングでアクセスすることによって、競合が発生し、データが不整合になることを防ぐために使います。

使用上の注意点

ブロック内の全スレッドが同時に呼び出す必要がある： __syncthreads()はブロック内の全てのスレッドが同時に呼び出す必要があります。もし、一部のスレッドが__syncthreads()に到達しない場合、デッドロックが発生し、プログラムが停止してしまいます。

ブロック外のスレッド間での同期はできない： __syncthreads()は、同じブロック内でのみ有効です。異なるブロック間での同期を行うためには、別の方法を用いる必要があります。

まとめ

__syncthreads()は、ブロック内のスレッドが協調して作業を行う際に非常に重要な関数です。特に、共有メモリを使用してスレッド間でデータをやり取りする際に、この同期機能を使うことで、正しいデータの整合性を保ちながら並列処理を効率的に進めることができます。

次は、スレッド間でデータを共有するために使う共有メモリ（shared）について、より詳しく見ていきます。

10. 共有メモリ shared の使い方と利点

CUDAプログラムで高いパフォーマンスを実現するために、共有メモリ（shared）は非常に重要な役割を果たします。共有メモリは、同じブロック内の全てのスレッドが共同で利用できる高速なメモリ領域です。CUDAのグローバルメモリに比べて非常に高速にアクセスできるため、共有メモリを効果的に活用することで、データの読み書き速度を大幅に向上させることができます。

共有メモリの特徴

ブロック内のスレッドで共有： 共有メモリは、同じブロック内のスレッド間でのみ共有されます。異なるブロック間ではデータを共有することができません。

• 高速なアクセス： 共有メモリは、グローバルメモリやテクスチャメモリよりも高速にアクセスできます。頻繁にアクセスするデータを共有メモリに保持することで、パフォーマンスを向上させることができます。
• サイズの制限： 共有メモリの使用量には制限があります。通常、1ブロックあたり数KB程度の共有メモリしか使えません。そのため、大規模なデータを処理する際は、共有メモリの使用量を工夫する必要があります。

shared の基本的な使い方

共有メモリは、カーネル内で__shared__キーワードを使って宣言します。次のコードは、共有メモリを使って配列の各要素に定数を加算する簡単な例です。

import pycuda.autoinit
import pycuda.driver as cuda
import numpy as np
from pycuda.compiler import SourceModule

# CUDAカーネルの定義
mod = SourceModule("""
__global__ void add_to_shared_memory(float *input, float *output)
{
    __shared__ float shared_data[10];  // 共有メモリの宣言
    int idx = threadIdx.x;

    // グローバルメモリから共有メモリへデータをコピー
    shared_data[idx] = input[idx];

    // 全スレッドを同期
    __syncthreads();

    // 共有メモリのデータを使用して計算
    output[idx] = shared_data[idx] + 10.0f;
}
""")

# カーネル関数を取得
add_to_shared_memory = mod.get_function("add_to_shared_memory")

# NumPy配列を用意
input_data = np.random.randn(10).astype(np.float32)

# GPUメモリに転送
input_gpu = cuda.mem_alloc(input_data.nbytes)
output_gpu = cuda.mem_alloc(input_data.nbytes)
cuda.memcpy_htod(input_gpu, input_data)

# カーネルを実行
add_to_shared_memory(input_gpu, output_gpu, block=(10, 1, 1))

# 結果を取得して表示
output_data = np.empty_like(input_data)
cuda.memcpy_dtoh(output_data, output_gpu)

print("入力データ:", input_data)
print("結果データ:", output_data)

例の解説

共有メモリの宣言： __shared__キーワードを使って、shared_dataという名前の共有メモリを宣言しています。この共有メモリは、各スレッドが共同で利用します。この例では、サイズが10の配列を共有メモリとして定義しています。

• 共有メモリへのデータコピー： 各スレッドは、グローバルメモリ（input）から自分の担当するデータを共有メモリにコピーしています。
• スレッドの同期： 全てのスレッドがデータの書き込みを終えるまで待機するために、__syncthreads()を使っています。この同期により、全スレッドが正しいデータにアクセスできる状態が保証されます。
• 共有メモリを使った計算： 同期が完了した後、各スレッドは共有メモリのデータを使用して計算を行い、その結果をグローバルメモリ（output）に書き戻しています。

共有メモリの応用例：畳み込みフィルタ

画像処理などの畳み込み演算においても、共有メモリは非常に有用です。畳み込みフィルタでは、各スレッドが自分の担当ピクセルとその周辺のピクセルにアクセスするため、同じデータを何度も読み込む必要があります。この際、共有メモリを使用してデータを一度読み込み、スレッド間で共有することでメモリアクセスのオーバーヘッドを削減できます。

以下に、3x3の畳み込みフィルタを実装する例を示します。

# 畳み込みカーネルの定義
mod = SourceModule("""
__global__ void convolve(float *image, float *output, int width, int height)
{
    __shared__ float shared_image[18][18];  // 16x16のブロックに余白を持たせる

    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;

    if (x >= width || y >= height) return;

    // グローバルメモリから共有メモリへデータをコピー
    shared_image[threadIdx.y + 1][threadIdx.x + 1] = image[y * width + x];

    // 共有メモリの端の処理
    if (threadIdx.x == 0 && x > 0)
        shared_image[threadIdx.y + 1][0] = image[y * width + (x - 1)];
    if (threadIdx.x == blockDim.x - 1 && x < width - 1)
        shared_image[threadIdx.y + 1][blockDim.x + 1] = image[y * width + (x + 1)];

    __syncthreads();

    // 畳み込みフィルタの適用
    float result = 0.0;
    result += shared_image[threadIdx.y][threadIdx.x];
    result += shared_image[threadIdx.y + 1][threadIdx.x];
    result += shared_image[threadIdx.y + 2][threadIdx.x];
    // 残りのフィルタ処理...

    output[y * width + x] = result;
}

共有メモリを使うメリット

• 高速化：グローバルメモリに直接アクセスするよりもはるかに高速な共有メモリを使うことで、スレッド間でデータを共有しながら高速な演算が可能になります。
• メモリアクセスの削減：複数のスレッドが同じデータを何度も読み取る場合、共有メモリに一度データを読み込んでおけば、その後は同じデータに高速にアクセスできます。

まとめ

共有メモリ（shared）は、同じブロック内のスレッドが共同で使える非常に高速なメモリです。効率的なデータ共有を可能にし、メモリアクセスを減らすことで、並列処理のパフォーマンスを大幅に向上させることができます。ただし、サイズに制約があるため、どのデータを共有メモリに格納するかを慎重に設計する必要があります。

次は、CUDAのストリーム機能を使って非同期処理を効率化する方法について紹介します。

11. CUDAストリームの基本と応用

CUDAストリーム（CUDA Streams）は、CUDAで非同期にタスクを実行するための重要な機能です。通常、カーネルの実行やメモリの転送は順番に行われますが、ストリームを使用することで、非同期に複数のタスクを同時並行で実行することが可能です。これにより、GPUのリソースを効率的に使い、パフォーマンスの向上が期待できます。

ストリームとは？

CUDAストリームは、一連の操作（カーネルの実行やメモリ転送）を順次実行するタスクキューのようなものです。CUDAでは、デフォルトで1つのストリームしか使用されませんが、複数のストリームを作成することで、メモリ転送とカーネル実行を非同期に行うことができるようになります。これにより、以下のようなメリットが得られます：

• 非同期実行：メモリの転送とカーネルの実行を同時に行う。
• 並列処理：複数のカーネルを異なるストリームで同時に実行する。

ストリームの使い方

• ストリームの作成： ストリームはcuda.Stream()で作成します。作成されたストリームは、カーネル実行やメモリ転送時に渡すことで、そのタスクを非同期で実行できます。
• ストリームにタスクを割り当てる： カーネルやメモリ転送を実行する際に、stream引数で使用するストリームを指定します。これにより、指定されたストリーム内でその操作が非同期に実行されます。
• ストリームの同期： cuda.Stream.synchronize()を使って、ストリーム内のすべての操作が完了するのを待機することができます。これにより、必要なタイミングで結果が揃っていることが保証されます。

例：ストリームを使った非同期カーネル実行

次に、ストリームを使ってメモリ転送とカーネルの実行を非同期で行う例を見てみましょう。

import pycuda.autoinit
import pycuda.driver as cuda
import pycuda.gpuarray as gpuarray
import numpy as np
from pycuda.compiler import SourceModule

# CUDAカーネルの定義
mod = SourceModule("""
__global__ void add_kernel(float *a, float *b, float *c, int n)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n)
        c[idx] = a[idx] + b[idx];
}
""")

# カーネル関数を取得
add_kernel = mod.get_function("add_kernel")

# データを準備
N = 1000000
a = np.random.randn(N).astype(np.float32)
b = np.random.randn(N).astype(np.float32)
c = np.empty_like(a)

# GPUメモリに転送
a_gpu = cuda.mem_alloc(a.nbytes)
b_gpu = cuda.mem_alloc(b.nbytes)
c_gpu = cuda.mem_alloc(c.nbytes)

# ストリームを作成
stream1 = cuda.Stream()
stream2 = cuda.Stream()

# 非同期にメモリ転送（2つのストリームを使う）
cuda.memcpy_htod_async(a_gpu, a, stream=stream1)
cuda.memcpy_htod_async(b_gpu, b, stream=stream2)

# 非同期にカーネルを実行（stream1を使用）
block_size = 256
grid_size = (N + block_size - 1) // block_size
add_kernel(a_gpu, b_gpu, c_gpu, np.int32(N), block=(block_size, 1, 1), grid=(grid_size, 1, 1), stream=stream1)

# 非同期にメモリ転送（結果をホストに戻す）
cuda.memcpy_dtoh_async(c, c_gpu, stream=stream1)

# ストリームの同期（stream1が終了するのを待つ）
stream1.synchronize()

# 結果を表示
print("計算結果の一部:", c[:10])

# ストリームを破棄
stream1 = None
stream2 = None

この例のポイント

• 複数のストリームの使用： ストリームを使うことで、メモリ転送とカーネル実行を非同期に行っています。stream1とstream2という2つのストリームを作成し、それぞれにタスクを割り当てています。
• 非同期メモリ転送： cuda.memcpy_htod_async()やcuda.memcpy_dtoh_async()は、非同期にメモリの転送を行います。これにより、CPU側で待機することなく、データ転送とカーネル実行を並行して行えます。
• ストリームの同期： stream1.synchronize()を使って、stream1の中で行われた全ての操作が完了するまで待機しています。この操作により、メモリ転送やカーネル実行が完全に終わってから次の処理に進むことが保証されます。

CUDAストリームの活用場面

• カーネルとメモリ転送の重ね合わせ： カーネルの実行とメモリ転送を非同期に実行することで、両方の操作を並列化できます。これにより、GPUの演算リソースとメモリバスの両方を効率的に利用することが可能になります。
• 複数のカーネルの同時実行： ストリームを使って複数のカーネルを並行して実行することで、GPUのマルチスレッディング機能を最大限に活用できます。これにより、特に複数の異なるタスクを同時に処理する場合に有効です。

まとめ

CUDAストリームを活用することで、非同期にカーネルを実行したり、メモリ転送を並列に行ったりすることが可能になります。これにより、CUDAプログラムのパフォーマンスを大幅に向上させることができます。ストリームを使うことで、タスクを重ね合わせたり、GPUリソースを最大限に活用できるため、大規模データ処理やリアルタイム処理において強力なツールとなります。