Nvidia製GPUのCudaを使ってジュリア集合を描く

Posted at 2024-01-21

Nvidia製GPUのCudaをつかうとCPUには到底不可能なほどの並列処理が可能となる。
当然計算処理が早くなるわけだが、両手をあげて喜んでいるわけには行かない。
CPUが利用するメモリとGPUが利用するメモリは切り離されていて相互にデーター転送が必要になる。
この転送速度がCuda処理よりもボトルネックになる場合はCPUで処理したほうが早いことになる。
Cudaを使うとなんでもかんでも処理速度が向上するわけでなないのだ。
サンプル例として取り上げられることが多い、２個のベクトルの足し算などはCudaで処理するほうが遅い場合が多い。
一方でメモリ転送量が少なく、計算の比重が大きい場合GPUで処理する効果が大きい。
今回はそんな例としてジュリア集合をCPUとCudaを使って速度比較してみた。

LinuxでのCuda開発方法　（Ubuntu)

Cuda(GPU)側のコンパイルはg++などではできないためNvidiaが用意したNVCCでコンパイルが必要である。

必要なものは以下

NVIDIA Driver
CUDA Toolkit

インストール方法などは他のサイトをみていただくほうがより良いと思われるのでここには記述しない。

CPUで処理する場合とCuda(GPU)の違い

ジュリア集合は複素平面（あとで複素平面がPNG画像になるのだが）の各複素数を漸化式を使って計算し、その値が収束するかしないかでその複素数が存在する位置の色を決めて画像を生成する。
この場合複素平面を表す画像のメモリ転送速度よりも各点における漸化式収束計算の方が時間がかかると思われるのでGPUで処理したほうが処理速度が早くなるはずである。

複素平面と漸化式の関係は以前記述した以下の記事を参考してしていただくとよい（手前味噌）
https://qiita.com/Guarneri001/items/1129e74e62017a792352

CPUとCuda(GPU)のコードの違い

早速だが、CPUとCuda(GPU)のコードの違いを見ていただきたい。

//CPU
#include <new>
#include <complex>
#include <cstddef>
#include <iostream>

inline float fminf(float a, float b)
{
    return a < b ? a : b;
}

inline float fmaxf(float a, float b)
{
    return a > b ? a : b;
}

inline float clamp(float f, float a, float b)
{
    return fmaxf(a, fminf(f, b));
}

// https://scrapbox.io/ePi5131/%E6%8B%A1%E5%BC%B5%E7%B7%A8%E9%9B%86%E3%81%AEHSV%E9%96%A2%E6%95%B0%E3%81%AEHSV%3ERGB%E5%A4%89%E6%8F%9B%E3%81%AE%E5%86%85%E5%AE%B9
void hsv2rgb_cpu(int h, int s, int v, int *r, int *g, int *b)
{
    if (h < 0)
        h += (1 - h / 360) * 360;
    if (360 < h)
        h %= 360;
    auto h1 = (h * 4096 + 50) / 120;
    auto s1 = (s * 4096 + 50) / 100;
    auto v1 = (v * 4096 + 50) / 100;
    auto h2 = h1 % 4096;
    auto a1{0}, a2{0};
    if (h2 < 2048)
    {
        a1 = (4096 - (2048 - h2) * s1 / 2048) * v1 / 4096;
        a2 = v1;
    }
    else
    {
        a2 = (4096 - (h2 - 2048) * s1 / 2048) * v1 / 4096;
        a1 = v1;
    }

    auto b1 = clamp((a2 * 255 + 2048) / 4096, 0, 255);
    auto b2 = clamp((a1 * 255 + 2048) / 4096, 0, 255);
    auto b3 = clamp(((4096 - s1) * v1 / 4096 * 255 + 2048) / 4096, 0, 255);

    switch (h1 / 4096)
    {
    case 1:
        *g = b1;
        *b = b2;
        *r = b3;
        break;
    case 2:
        *b = b1;
        *r = b2;
        *g = b3;
        break;
    default:
        *r = b1;
        *g = b2;
        *b = b3;
        break;
    }
}

int julia_CPU(int x, int y, int view_size)
{
    const auto scale = 1.5f;
    auto jx = scale * (float)(view_size / 2 - x) / (view_size / 2);
    auto jy = scale * (float)(view_size / 2 - y) / (view_size / 2);

    std::complex<float> c(-0.8f, 0.156f);
    std::complex<float> z(jx, jy);

    auto i = 0;
    for (i = 0; i < 360; i++)
    {
        z = z * z + c;
        if (std::norm(z) > 1000)
            break;
    }
    return i;
}

extern "C" unsigned char *JuliaCPU(std::size_t size, int view_size)
{
    auto *ptr = new unsigned char[size];
    auto r{0}, g{0}, b{0};

    for (auto y = 0; y < view_size; y++)
    {
        for (auto x = 0; x < view_size; x++)
        {
            auto offset = x + y * view_size;
            auto value = julia_CPU(x, y, view_size);

            if (value >= 0 && value <= 20)
            {
                ptr[offset * 4 + 0] = 255;
                ptr[offset * 4 + 1] = 255;
                ptr[offset * 4 + 2] = 255;
                ptr[offset * 4 + 3] = 255;
                continue;
            }

            hsv2rgb_cpu(value, 100, 100, &r, &g, &b);

            ptr[offset * 4 + 0] = r;
            ptr[offset * 4 + 1] = g;
            ptr[offset * 4 + 2] = b;
            ptr[offset * 4 + 3] = 255;
        }
    }
    return ptr;
}

CPUコードでは x軸とy軸の２重ループになっている。

//GPU
#include <iostream>
#include <cstddef>
#include <thrust/complex.h>
#include <cuda.h>

constexpr auto threads_perblock = 1024;

static void HandleError(cudaError_t err, const char *file, int line)
{
    if (err != cudaSuccess)
    {
        std::cout << cudaGetErrorString(err) << "  " << file << "  " << line << std::endl;
        exit(EXIT_FAILURE);
    }
}
#define HANDLE_ERROR(err) (HandleError(err, __FILE__, __LINE__))

inline __device__ float fminf(float a, float b)
{
    return a < b ? a : b;
}

inline __device__ float fmaxf(float a, float b)
{
    return a > b ? a : b;
}

inline __device__ float clamp(float f, float a, float b)
{
    return fmaxf(a, fminf(f, b));
}

__device__ void hsv2rgb_gpu(int h, int s, int v, int *r, int *g, int *b)
{
    if (h < 0)
        h += (1 - h / 360) * 360;
    if (360 < h)
        h %= 360;
    auto h1 = (h * 4096 + 50) / 120;
    auto s1 = (s * 4096 + 50) / 100;
    auto v1 = (v * 4096 + 50) / 100;
    auto h2 = h1 % 4096;
    auto a1{0}, a2{0};
    if (h2 < 2048)
    {
        a1 = (4096 - (2048 - h2) * s1 / 2048) * v1 / 4096;
        a2 = v1;
    }
    else
    {
        a2 = (4096 - (h2 - 2048) * s1 / 2048) * v1 / 4096;
        a1 = v1;
    }

    auto b1 = clamp((a2 * 255 + 2048) / 4096, 0, 255);
    auto b2 = clamp((a1 * 255 + 2048) / 4096, 0, 255);
    auto b3 = clamp(((4096 - s1) * v1 / 4096 * 255 + 2048) / 4096, 0, 255);

    switch (h1 / 4096)
    {
    case 1:
        *g = b1;
        *b = b2;
        *r = b3;
        break;
    case 2:
        *b = b1;
        *r = b2;
        *g = b3;
        break;
    default:
        *r = b1;
        *g = b2;
        *b = b3;
        break;
    }
}

__device__ int julia_cuda(int x, int y, int view_size)
{
    const auto scale = 1.5f;
    auto jx = scale * static_cast<float>(view_size / 2 - x) / (view_size / 2);
    auto jy = scale * static_cast<float>(view_size / 2 - y) / (view_size / 2);

    thrust::complex<float> c(-0.8f, 0.156f);
    thrust::complex<float> z(jx, jy);

    auto i = 0;
    for (i = 0; i < 360; i++)
    {
        z = z * z + c;
        if (thrust::norm(z) > 1000)
            break;
    }

    return i;
}

__global__ void kernel(unsigned char *ptr, int view_size)
{
    auto r{0}, g{0}, b{0};
    auto x = blockIdx.x;
    auto y = threadIdx.x;
    auto offset = x + y * gridDim.x;
    auto value = julia_cuda(x, y, view_size);

    hsv2rgb_gpu(value, 100, 100, &r, &g, &b);

    if (value >= 0 && value <= 20)
    {
        ptr[offset * 4 + 0] = 255;
        ptr[offset * 4 + 1] = 255;
        ptr[offset * 4 + 2] = 255;
        ptr[offset * 4 + 3] = 255;
    }
    else
    {
        ptr[offset * 4 + 0] = r;
        ptr[offset * 4 + 1] = g;
        ptr[offset * 4 + 2] = b;
        ptr[offset * 4 + 3] = 255;
    }
}

extern "C" unsigned char *JuliaGPU(std::size_t size, int view_size)
{
    unsigned char *ptr_gpu;

    HANDLE_ERROR(cudaMalloc((void **)&ptr_gpu, size));

    auto *ptr = new unsigned char[size];
    HANDLE_ERROR(cudaMemcpy(ptr_gpu, ptr, size, cudaMemcpyHostToDevice));

    auto blocks_per_grid = ((view_size * view_size) + threads_perblock - 1) / threads_perblock;
    std::cout << "CUDA kernel [" << blocks_per_grid << "] blocks [" << threads_perblock << "] threads" << std::endl;
    kernel<<<blocks_per_grid, threads_perblock>>>(ptr_gpu, view_size);

    HANDLE_ERROR(cudaMemcpy(ptr, ptr_gpu, size, cudaMemcpyDeviceToHost));
    HANDLE_ERROR(cudaFree(ptr_gpu));
    return ptr;
}

Cuda(GPU)コードでは x軸とy軸の２重ループが無くなり kernel関数が追加され、CPU側からのメモリ転送用コードが記述されている。CPUでループ処理していた部分はCuda(GPU)コードではカーネル処理に変更され、このカーネルコードを各Cudaコアが並列に処理することによりCPUを上回る並列数で処理される。

上記の見慣れない記述 <<<blocks_per_grid, threads_perblock>>> があるが、これがCudaにおける並列処理用のパラメータでグリッド内のブロック数とそのブロック内で実行されるスレッド数を表している。

今回は　blocks_per_grid = 1024, threads_perblock = 1024 で設定、PNG画像を　1024x1024　で試してみている blockIdx.x が x 軸、threadIdx.x が y軸としてカーネルに渡され各担当の位置の計算をCudaコアで行っている。

処理速度結果

試したPCスペックは以下

CPU i5-4590
GPU GTX970

CPU total result.. 1755.27 ms
GPU total result.. 80.3649 ms

Cudaの方が20倍も早かった。

両方とも全く同じ画像が出力された。（当たり前か）

以下に完全なソースコードを公開
https://github.com/Guarneri009/cuda-julia

CUDA_ARCHITECTURES "52" の部分は搭載されているGPUに合わせて調整が必要かもしれない

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

Nvidia製GPUのCudaを使ってジュリア集合を描く

LinuxでのCuda開発方法 （Ubuntu)

CPUで処理する場合とCuda(GPU)の違い

CPUとCuda(GPU)のコードの違い

処理速度結果

LinuxでのCuda開発方法　（Ubuntu)