PyCUDAで始めるGPUプログラミング入門

第1章: PyCUDAとは何か

PyCUDAは、PythonからNVIDIA GPUの計算能力を活用するためのツールキットです。CUDAプログラミングの知識をPythonの使いやすさと組み合わせることで、高性能な並列計算を実現できます。

PyCUDAを使うと、以下のようなメリットがあります：

1. Pythonの柔軟性とGPUの高速処理を組み合わせられる
2. メモリ管理が自動化され、エラーが起こりにくい
3. CUDAカーネルを動的に生成・コンパイルできる

では、最初の簡単なPyCUDAプログラムを見てみましょう：

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

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

# データの準備
a = np.random.randn(400).astype(np.float32)
a_gpu = cuda.mem_alloc(a.nbytes)
cuda.memcpy_htod(a_gpu, a)

# カーネルの実行
func = mod.get_function("multiply_by_two")
func(a_gpu, block=(400,1,1), grid=(1,1))

# 結果の取得
a_doubled = np.empty_like(a)
cuda.memcpy_dtoh(a_doubled, a_gpu)

print("元の配列:", a)
print("2倍にした配列:", a_doubled)

このコードは、400要素の浮動小数点数配列を生成し、それぞれの要素を2倍にするGPU関数を実行しています。

第2章: PyCUDAのセットアップ

PyCUDAを使用するには、まずNVIDIA GPUとCUDAツールキットがインストールされている必要があります。その後、以下の手順でPyCUDAをインストールします：

1. Anacondaを使用している場合：

   conda install -c conda-forge pycuda
   

2. pipを使用する場合：

   pip install pycuda
   

インストールが完了したら、以下のコードでPyCUDAが正しく動作するか確認しましょう：

import pycuda.autoinit
import pycuda.driver as cuda

print(f"CUDA version: {cuda.get_version()}")
print(f"Device name: {cuda.Device(0).name()}")

このコードはCUDAのバージョンとデバイス名を表示します。

第3章: 基本的なベクトル演算

GPUを使って基本的なベクトル演算を行ってみましょう。ここでは、二つのベクトルの要素ごとの加算を行います。

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

mod = SourceModule("""
__global__ void add_vectors(float *a, float *b, float *result, int n)
{
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if (idx < n) {
        result[idx] = a[idx] + b[idx];
    }
}
""")

def vector_add(a, b):
    n = a.shape[0]
    # GPUメモリの確保
    a_gpu = cuda.mem_alloc(a.nbytes)
    b_gpu = cuda.mem_alloc(b.nbytes)
    result_gpu = cuda.mem_alloc(a.nbytes)

    # データのGPUへの転送
    cuda.memcpy_htod(a_gpu, a)
    cuda.memcpy_htod(b_gpu, b)

    # カーネルの実行
    func = mod.get_function("add_vectors")
    block_size = 256
    grid_size = (n + block_size - 1) // block_size
    func(a_gpu, b_gpu, result_gpu, np.int32(n),
         block=(block_size, 1, 1), grid=(grid_size, 1))

    # 結果の取得
    result = np.empty_like(a)
    cuda.memcpy_dtoh(result, result_gpu)

    return result

# テスト
a = np.random.randn(1000000).astype(np.float32)
b = np.random.randn(1000000).astype(np.float32)

result = vector_add(a, b)
print("GPU結果の最初の10要素:", result[:10])
print("CPU結果の最初の10要素:", (a + b)[:10])
print("結果が一致:", np.allclose(result, a + b))

このコードでは、100万要素のベクトルの加算をGPUで行っています。CPUでの計算結果と比較して、正しく計算できていることを確認しています。

第4章: 行列乗算

行列乗算は、多くの科学技術計算で重要な演算です。GPUを使うことで、大規模な行列乗算を高速に処理できます。

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

mod = SourceModule("""
__global__ void matrix_multiply(float *a, float *b, float *c, int m, int n, int k)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    float sum = 0.0f;
    if (row < m && col < k) {
        for (int i = 0; i < n; i++) {
            sum += a[row * n + i] * b[i * k + col];
        }
        c[row * k + col] = sum;
    }
}
""")

def gpu_matrix_multiply(a, b):
    m, n = a.shape
    n, k = b.shape
    
    # GPUメモリの確保
    a_gpu = cuda.mem_alloc(a.nbytes)
    b_gpu = cuda.mem_alloc(b.nbytes)
    c_gpu = cuda.mem_alloc(m * k * 4)  # float32なので4バイト

    # データのGPUへの転送
    cuda.memcpy_htod(a_gpu, a)
    cuda.memcpy_htod(b_gpu, b)

    # カーネルの実行
    func = mod.get_function("matrix_multiply")
    block = (16, 16, 1)
    grid = ((k + block[0] - 1) // block[0], (m + block[1] - 1) // block[1])
    func(a_gpu, b_gpu, c_gpu, np.int32(m), np.int32(n), np.int32(k),
         block=block, grid=grid)

    # 結果の取得
    c = np.empty((m, k), dtype=np.float32)
    cuda.memcpy_dtoh(c, c_gpu)

    return c

# テスト
a = np.random.randn(1000, 1000).astype(np.float32)
b = np.random.randn(1000, 1000).astype(np.float32)

c_gpu = gpu_matrix_multiply(a, b)
c_cpu = np.dot(a, b)

print("GPU結果の一部:")
print(c_gpu[:5, :5])
print("\nCPU結果の一部:")
print(c_cpu[:5, :5])
print("\n結果が一致:", np.allclose(c_gpu, c_cpu, rtol=1e-5))

このコードでは、1000x1000の行列同士の乗算をGPUで行っています。大規模な行列でも高速に計算できることがわかります。

第5章: 画像処理：ぼかしフィルタ

GPUは画像処理にも非常に適しています。ここでは、簡単なぼかしフィルタを実装してみましょう。

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

mod = SourceModule("""
__global__ void blur_filter(unsigned char *in, unsigned char *out, int width, int height)
{
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (col < width && row < height) {
        int pixVal = 0;
        int pixels = 0;
        for (int blurrow = -2; blurrow < 3; blurrow++) {
            for (int blurcol = -2; blurcol < 3; blurcol++) {
                int currow = row + blurrow;
                int curcol = col + blurcol;
                if (currow > -1 && currow < height && curcol > -1 && curcol < width) {
                    pixVal += in[currow * width + curcol];
                    pixels++;
                }
            }
        }
        out[row * width + col] = (unsigned char)(pixVal / pixels);
    }
}
""")

def gpu_blur(image):
    img_array = np.array(image, dtype=np.uint8)
    height, width = img_array.shape

    # GPUメモリの確保
    in_gpu = cuda.mem_alloc(img_array.nbytes)
    out_gpu = cuda.mem_alloc(img_array.nbytes)

    # データのGPUへの転送
    cuda.memcpy_htod(in_gpu, img_array)

    # カーネルの実行
    func = mod.get_function("blur_filter")
    block = (16, 16, 1)
    grid = ((width + block[0] - 1) // block[0], (height + block[1] - 1) // block[1])
    func(in_gpu, out_gpu, np.int32(width), np.int32(height),
         block=block, grid=grid)

    # 結果の取得
    result = np.empty_like(img_array)
    cuda.memcpy_dtoh(result, out_gpu)

    return Image.fromarray(result)

# テスト
image = Image.open("input_image.jpg").convert("L")  # グレースケールに変換
blurred = gpu_blur(image)
blurred.save("output_blurred.jpg")

print("ぼかし処理が完了しました。output_blurred.jpgを確認してください。")

このコードでは、入力画像に5x5のぼかしフィルタを適用しています。大きな画像でも高速に処理できるのがGPUの強みです。

第6章: 並列リダクション：最大値の計算

並列リダクションは、大規模なデータセットから単一の結果を得るための重要な技術です。ここでは、配列の最大値を求める並列リダクションを実装します。

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

mod = SourceModule("""
__global__ void reduce_max(float *in, float *out, int N)
{
    __shared__ float sdata[256];
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    
    sdata[tid] = (i < N) ? in[i] : -INFINITY;
    __syncthreads();

    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] = max(sdata[tid], sdata[tid + s]);
        }
        __syncthreads();
    }

    if (tid == 0) out[blockIdx.x] = sdata[0];
}
""")

def gpu_max(arr):
    N = arr.shape[0]
    
    # GPUメモリの確保
    in_gpu = cuda.mem_alloc(arr.nbytes)
    cuda.memcpy_htod(in_gpu, arr)

    # ブロックサイズとグリッドサイズの設定
    block_size = 256
    grid_size = (N + block_size - 1) // block_size
    out_gpu = cuda.mem_alloc(grid_size * 4)  # float32なので4バイト

    # カーネルの実行
    func = mod.get_function("reduce_max")
    func(in_gpu, out_gpu, np.int32(N), block=(block_size, 1, 1), grid=(grid_size, 1))

    # 結果の取得
    result = np.empty(grid_size, dtype=np.float32)
    cuda.memcpy_dtoh(result, out_gpu)

    # 最終的な最大値を計算
    return np.max(result)

# テスト
arr = np.random.randn(1000000).astype(np.float32)
gpu_result = gpu_max(arr)
cpu_result = np.max(arr)

print(f"GPU結果: {gpu_result}")
print(f"CPU結果: {cpu_result}")
print(f"結果が一致: {np.isclose(gpu_result, cpu_result)}")

このコードでは、100万要素の配列から最大値を求めています。GPUの並列処理能力を活かすことで、大規模なデータセットでも高速に計算できます。

第7章: カスタムデータ型の使用

PyCUDAでは、カスタムデータ型を定義して使用することができます。これにより、複雑なデータ構造を扱うことが可能になります。

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

mod = SourceModule("""
struct Complex {
    float real;
    float imag;
};

__global__ void add_complex(Complex *a, Complex *b, Complex *result, int n)
{
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if (idx < n) {
        result[idx].real = a[idx].real + b[idx].real;
        result[idx].imag = a[idx].imag + b[idx].imag;
    }
}
""")

def complex_add(a, b):
    n = len(a)
    
    # NumPy構造体配列の作成
    complex_dtype = np.dtype([('real', np.float32), ('imag', np.float32)])
    a_array = np.array([(c.real, c.imag) for c in a], dtype=complex_dtype)
    b_array = np.array([(c.real, c.imag) for c in b], dtype=complex_dtype)
    
    # GPUメモリの確保
    a_gpu = cuda.mem_alloc(a_array.nbytes)
    b_gpu = cuda.mem_alloc(b_array.nbytes)
    result_gpu = cuda.mem_alloc(a_array.nbytes)

    # データのGPUへの転送
    cuda.memcpy_htod(a_gpu, a_array)
    cuda.memcpy_htod(b_gpu, b_array)

    # カーネルの実行
    func = mod.get_function("add_complex")
    block_size = 256
    grid_size = (n + block_size - 1) // block_size
    func(a_gpu, b_gpu, result_gpu, np.int32(n),
         block=(block_size, 1, 1), grid=(grid_size, 1))

    # 結果の取得
    result_array = np.empty(n, dtype=complex_dtype)
    cuda.memcpy_dtoh(result_array, result_gpu)

    # 複素数オブジェクトに変換して返す
    return [complex(x['real'], x['imag']) for x in result_array]

# テスト
a = [complex(1, 2), complex(3, 4), complex(5, 6)]
b = [complex(7, 8), complex(9, 10), complex(11, 12)]

result = complex_add(a, b)
print("GPU結果:", result)
print("CPU結果:", [a[i] + b[i] for i in range(len(a))])

このコードでは、複素数のカスタムデータ型を定義し、複素数の配列同士の加算を行っています。カスタムデータ型を使用することで、より複雑なデータ構造をGPUで効率的に処理できます。

第8章: スレッド同期とシェアードメモリ

GPUプログラミングでは、スレッド間の同期とシェアードメモリの使用が重要です。ここでは、行列の転置を例に、これらの概念を説明します。

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

mod = SourceModule("""
#define BLOCK_SIZE 16

__global__ void transpose(float *in, float *out, int width, int height)
{
    __shared__ float tile[BLOCK_SIZE][BLOCK_SIZE+1];
    
    int x = blockIdx.x * BLOCK_SIZE + threadIdx.x;
    int y = blockIdx.y * BLOCK_SIZE + threadIdx.y;
    
    if (x < width && y < height) {
        tile[threadIdx.y][threadIdx.x] = in[y * width + x];
    }
    
    __syncthreads();
    
    x = blockIdx.y * BLOCK_SIZE + threadIdx.x;
    y = blockIdx.x * BLOCK_SIZE + threadIdx.y;
    
    if (x < height && y < width) {
        out[y * height + x] = tile[threadIdx.x][threadIdx.y];
    }
}
""")

def gpu_transpose(matrix):
    height, width = matrix.shape
    
    # GPUメモリの確保
    in_gpu = cuda.mem_alloc(matrix.nbytes)
    out_gpu = cuda.mem_alloc(matrix.nbytes)

    # データのGPUへの転送
    cuda.memcpy_htod(in_gpu, matrix)

    # カーネルの実行
    func = mod.get_function("transpose")
    block = (16, 16, 1)
    grid = ((width + block - 1) // block, (height + block - 1) // block)
    func(in_gpu, out_gpu, np.int32(width), np.int32(height),
         block=block, grid=grid)

    # 結果の取得
    result = np.empty((width, height), dtype=np.float32)
    cuda.memcpy_dtoh(result, out_gpu)

    return result

# テスト
matrix = np.random.randn(1024, 1024).astype(np.float32)
transposed = gpu_transpose(matrix)

print("元の行列の形状:", matrix.shape)
print("転置後の行列の形状:", transposed.shape)
print("結果が正しいか:", np.allclose(transposed, matrix.T))

このコードでは、行列の転置を効率的に行うために、シェアードメモリとスレッド同期を使用しています。__syncthreads()を使ってブロック内のすべてのスレッドを同期させ、シェアードメモリを使って効率的にデータを転送しています。

第9章: ストリームを使った非同期処理

CUDAストリームを使用すると、複数の処理を非同期的に実行できます。これにより、GPUの利用効率を高めることができます。

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

mod = SourceModule("""
__global__ void vector_add(float *a, float *b, float *c, int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        c[i] = a[i] + b[i];
    }
}
""")

def async_vector_add(a, b, num_streams=4):
    n = len(a)
    chunk_size = n // num_streams

    # GPUメモリの確保
    a_gpu = cuda.mem_alloc(a.nbytes)
    b_gpu = cuda.mem_alloc(b.nbytes)
    c_gpu = cuda.mem_alloc(a.nbytes)

    # カーネル関数の取得
    func = mod.get_function("vector_add")

    # ストリームの作成
    streams = [cuda.Stream() for _ in range(num_streams)]

    for i in range(num_streams):
        start = i * chunk_size
        end = (i + 1) * chunk_size if i < num_streams - 1 else n

        # 非同期データ転送
        cuda.memcpy_htod_async(int(a_gpu) + start * 4, a[start:end], streams[i])
        cuda.memcpy_htod_async(int(b_gpu) + start * 4, b[start:end], streams[i])

        # 非同期カーネル実行
        func(a_gpu, b_gpu, c_gpu, np.int32(n),
             block=(256, 1, 1), grid=((end - start + 255) // 256, 1),
             stream=streams[i])

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

    return c

# テスト
n = 10000000
a = np.random.randn(n).astype(np.float32)
b = np.random.randn(n).astype(np.float32)

c = async_vector_add(a, b)

print("結果の一部:", c[:10])
print("正しく計算できているか:", np.allclose(c, a + b))

このコードでは、大きな配列を複数のストリームに分割して非同期的に処理しています。これにより、データ転送と計算のオーバーラップが可能になり、全体的な処理時間を短縮できます。

第10章: テクスチャメモリの活用

テクスチャメモリは、2次元データの読み取りに最適化されたGPUのメモリです。画像処理などで効果を発揮します。

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

mod = SourceModule("""
texture<float, 2, cudaReadModeElementType> tex;

__global__ void rotate_image(float *output, int width, int height, float angle)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    float cx = width / 2.0f;
    float cy = height / 2.0f;
    
    if (x < width && y < height) {
        float nx = cosf(angle) * (x - cx) - sinf(angle) * (y - cy) + cx;
        float ny = sinf(angle) * (x - cx) + cosf(angle) * (y - cy) + cy;
        
        if (nx >= 0 && nx < width && ny >= 0 && ny < height) {
            output[y * width + x] = tex2D(tex, nx, ny);
        } else {
            output[y * width + x] = 0.0f;
        }
    }
}
""")

def rotate_image_gpu(image, angle):
    # 画像をNumPy配列に変換
    img_array = np.array(image, dtype=np.float32) / 255.0
    height, width = img_array.shape

    # テクスチャの設定
    texref = mod.get_texref("tex")
    cuda.matrix_to_texref(img_array, texref, order="C")
    texref.set_flags(cuda.TRSF_NORMALIZED_COORDINATES)
    texref.set_filter_mode(cuda.filter_mode.LINEAR)
    texref.set_address_mode(0, cuda.address_mode.WRAP)
    texref.set_address_mode(1, cuda.address_mode.WRAP)

    # 出力用GPUメモリの確保
    output_gpu = cuda.mem_alloc(img_array.nbytes)

    # カーネルの実行
    func = mod.get_function("rotate_image")
    block = (16, 16, 1)
    grid = ((width + block - 1) // block, (height + block - 1) // block)
    func(output_gpu, np.int32(width), np.int32(height), np.float32(angle),
         texrefs=[texref], block=block, grid=grid)

    # 結果の取得
    output = np.empty_like(img_array)
    cuda.memcpy_dtoh(output, output_gpu)

    # 画像に変換して返す
    return Image.fromarray((output * 255).astype(np.uint8))

# テスト
image = Image.open("input_image.jpg").convert("L")  # グレースケールに変換
rotated = rotate_image_gpu(image, np.pi / 4)  # 45度回転
rotated.save("output_rotated.jpg")

print("画像の回転が完了しました。output_rotated.jpgを確認してください。")

このコードでは、テクスチャメモリを使用して画像を回転させています。テクスチャメモリを使用することで、2次元データの補間や境界処理が容易になります。

第11章: 動的並列処理

CUDA 5.0以降では、GPUカーネル内から別のカーネルを起動する動的並列処理が可能になりました。これにより、より複雑な並列アルゴリズムを実装できます。

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

mod = SourceModule("""
#include <stdio.h>

__global__ void child_kernel(int depth) {
    if (depth > 0) {
        printf("Depth %d: Thread %d in block %d\\n", depth, threadIdx.x, blockIdx.x);
        if (threadIdx.x == 0) {
            child_kernel<<<1, 32>>>(depth - 1);
        }
    }
}

__global__ void parent_kernel(int max_depth) {
    child_kernel<<<1, 32>>>(max_depth);
}
""", no_extern_c=True)

def run_dynamic_parallelism(max_depth):
    func = mod.get_function("parent_kernel")
    func(np.int32(max_depth), block=(1, 1, 1), grid=(1, 1))

# テスト
run_dynamic_parallelism(3)

このコードでは、親カーネルが子カーネルを呼び出し、さらに子カーネルが再帰的に自身を呼び出しています。これにより、GPUで木構造のようなアルゴリズムを実装できます。

第12章: アトミック操作

複数のスレッドが同時にメモリにアクセスする場合、アトミック操作が重要になります。ここでは、ヒストグラムの計算を例に説明します。

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

mod = SourceModule("""
__global__ void histogram(int *input, int *output, int n, int bins)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        int bin = input[idx] % bins;
        atomicAdd(&output[bin], 1);
    }
}
""")

def gpu_histogram(data, bins):
    n = len(data)
    
    # GPUメモリの確保
    input_gpu = cuda.mem_alloc(data.nbytes)
    output_gpu = cuda.mem_alloc(bins * 4)  # int32なので4バイト

    # データのGPUへの転送
    cuda.memcpy_htod(input_gpu, data)
    cuda.memcpy_htod(output_gpu, np.zeros(bins, dtype=np.int32))

    # カーネルの実行
    func = mod.get_function("histogram")
    block_size = 256
    grid_size = (n + block_size - 1) // block_size
    func(input_gpu, output_gpu, np.int32(n), np.int32(bins),
         block=(block_size, 1, 1), grid=(grid_size, 1))

    # 結果の取得
    result = np.empty(bins, dtype=np.int32)
    cuda.memcpy_dtoh(result, output_gpu)

    return result

# テスト
data = np.random.randint(0, 100, 1000000, dtype=np.int32)
bins = 100

gpu_result = gpu_histogram(data, bins)
cpu_result, _ = np.histogram(data, bins=bins, range=(0, bins))

print("GPU結果の一部:", gpu_result[:10])
print("CPU結果の一部:", cpu_result[:10])
print("結果が一致:", np.all(gpu_result == cpu_result))

このコードでは、atomicAddを使用してヒストグラムを計算しています。アトミック操作により、複数のスレッドが同時に同じビンをインクリメントしても、正確な結果が得られます。

第13章: 最適化テクニック

GPUプログラミングでは、パフォーマンスを最大化するためのさまざまな最適化テクニックがあります。ここでは、行列乗算を例に、いくつかの最適化テクニックを紹介します。

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

mod = SourceModule("""
#define TILE_WIDTH 32

__global__ void matrix_mul_optimized(float *A, float *B, float *C, int M, int N, int K)
{
    __shared__ float As[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Bs[TILE_WIDTH][TILE_WIDTH];
    
    int bx = blockIdx.x;  int by = blockIdx.y;
    int tx = threadIdx.x; int ty = threadIdx.y;
    
    int row = by * TILE_WIDTH + ty;
    int col = bx * TILE_WIDTH + tx;
    
    float sum = 0.0f;
    
    for (int t = 0; t < (N - 1) / TILE_WIDTH + 1; ++t) {
        if (row < M && t * TILE_WIDTH + tx < N)
            As[ty][tx] = A[row * N + t * TILE_WIDTH + tx];
        else
            As[ty][tx] = 0.0f;
        
        if (t * TILE_WIDTH + ty < N && col < K)
            Bs[ty][tx] = B[(t * TILE_WIDTH + ty) * K + col];
        else
            Bs[ty][tx] = 0.0f;
        
        __syncthreads();
        
        for (int i = 0; i < TILE_WIDTH; ++i)
            sum += As[ty][i] * Bs[i][tx];
        
        __syncthreads();
    }
    
    if (row < M && col < K)
        C[row * K + col] = sum;
}
""")

def gpu_matrix_multiply_optimized(A, B):
    M, N = A.shape
    N, K = B.shape
    
    # GPUメモリの確保
    A_gpu = cuda.mem_alloc(A.nbytes)
    B_gpu = cuda.mem_alloc(B.nbytes)
    C_gpu = cuda.mem_alloc(M * K * 4)  # float32なので4バイト

    # データのGPUへの転送
    cuda.memcpy_htod(A_gpu, A)
    cuda.memcpy_htod(B_gpu, B)

    # カーネルの実行
    func = mod.get_function("matrix_mul_optimized")
    tile_width = 32
    grid = ((K + tile_width - 1) // tile_width, (M + tile_width - 1) // tile_width)
    func(A_gpu, B_gpu, C_gpu, np.int32(M), np.int32(N), np.int32(K),
         block=(tile_width, tile_width, 1), grid=grid)

    # 結果の取得
    C = np.empty((M, K), dtype=np.float32)
    cuda.memcpy_dtoh(C, C_gpu)

    return C

# テスト
A = np.random.randn(1024, 1024).astype(np.float32)
B = np.random.randn(1024, 1024).astype(np.float32)

C_gpu = gpu_matrix_multiply_optimized(A, B)
C_cpu = np.dot(A, B)

print("GPU結果の一部:")
print(C_gpu[:5, :5])
print("\nCPU結果の一部:")
print(C_cpu[:5, :5])
print("\n結果が一致:", np.allclose(C_gpu, C_cpu, rtol=1e-5))

このコードでは、以下の最適化テクニックを使用しています：

1. シェアードメモリの使用：タイル状にデータを分割し、シェアードメモリにロードすることで、グローバルメモリへのアクセスを減らしています。
2. ループアンローリング：内部ループを展開することで、命令レベルの並列性を向上させています。
3. 境界チェック：タイルの境界でのチェックを行い、配列の範囲外アクセスを防いでいます。

これらの最適化により、単純な実装と比べて大幅なパフォーマンス向上が期待できます。

第14章: プロファイリングとデバッグ

GPUプログラムのパフォーマンスを最適化するには、プロファイリングが重要です。また、エラーを効率的に見つけるためのデバッグ技術も必要です。

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

mod = SourceModule("""
__global__ void vector_add(float *a, float *b, float *c, int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        c[i] = a[i] + b[i];
    }
}
""")

def profile_vector_add(n):
    a = np.random.randn(n).astype(np.float32)
    b = np.random.randn(n).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)

    # データのGPUへの転送
    cuda.memcpy_htod(a_gpu, a)
    cuda.memcpy_htod(b_gpu, b)

    # カーネルの実行
    func = mod.get_function("vector_add")
    
    # プロファイリングの開始
    start = cuda.Event()
    end = cuda.Event()
    start.record()

    func(a_gpu, b_gpu, c_gpu, np.int32(n),
         block=(256, 1, 1), grid=((n + 255) // 256, 1))

    # プロファイリングの終了
    end.record()
    end.synchronize()

    # 実行時間の計算
    time_ms = start.time_till(end)

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

    return time_ms, c

# テスト
n = 10000000
time_ms, result = profile_vector_add(n)

print(f"実行時間: {time_ms:.2f} ミリ秒")
print("結果の一部:", result[:10])

このコードでは、cuda.Eventを使用してカーネルの実行時間を測定しています。これにより、異なる実装や入力サイズでのパフォーマンスを比較できます。

デバッグについては、以下のテクニックが有効です：

1. printf文の使用：カーネル内でprintfを使用して、中間結果や変数の値を出力できます。
2. CUDA-GDB：NVIDIAが提供するデバッガを使用して、カーネルをステップ実行できます。
3. エラーチェック：cudaGetLastError()を使用して、カーネル実行後にエラーをチェックします。

第15章: 高度なトピック

最後に、PyCUDAを使用する上で知っておくと良い高度なトピックをいくつか紹介します。

1. マルチGPUプログラミング：
   複数のGPUを使用して並列処理を行うことができます。

import pycuda.driver as cuda

# 利用可能なGPUの数を取得
num_gpus = cuda.Device.count()
print(f"利用可能なGPU数: {num_gpus}")

# 各GPUの情報を表示
for i in range(num_gpus):
    gpu = cuda.Device(i)
    print(f"GPU {i}:")
    print(f"  名前: {gpu.name()}")
    print(f"  メモリ: {gpu.total_memory() // (1024**2)} MB")
    print(f"  コンピュート能力: {gpu.compute_capability()}")

1. 非同期メモリ転送：
   cuda.memcpy_htod_asyncとcuda.memcpy_dtoh_asyncを使用して、計算とメモリ転送をオーバーラップさせることができます。

2. 持続的カーネル：
   CUDA 8.0以降では、持続的カーネルを使用して、カーネルの起動オーバーヘッドを削減できます。

3. 動的並列処理：
   GPUカーネル内から別のカーネルを起動する動的並列処理を使用して、より複雑なアルゴリズムを実装できます。

4. 協調グループ：
   CUDA 9.0以降では、協調グループを使用して、より柔軟なスレッド同期とコミュニケーションが可能になりました。

これらの高度なトピックを理解し活用することで、より効率的で柔軟なGPUプログラミングが可能になります。

以上で、PyCUDAに関する15章構成の解説を終わります。この記事を通じて、PyCUDAの基本から応用まで幅広く学ぶことができたと思います。GPUプログラミングは奥が深く、常に新しい技術や最適化手法が登場しています。継続的に学習し、実践することで、より効率的なGPUプログラムを開発できるようになるでしょう。