超速の高速化エンジニアリング。CUDA GPU : 2.45秒 C++ コンパイラ最適化 : 2.99秒 AVX-512命令で最適化: 0.20秒。GPUをぶっちぎった。

タイトル: 超速の伝説

ある未来の世界で、計算の神秘を解き明かす科学者たちは、時空を超える速度を追い求めていた。その中に、「アレックス」という若き科学者がいた。彼は数式とコンピュータの力を駆使して、全宇宙の謎を解明する夢を抱いていた。

アレックスの最新のプロジェクトは、宇宙の始まりに隠された「マンデルブロ集合」の深奥なパターンを解明することだった。これは、無限の複雑さを持ち、未知の美しさを秘めた数式の集まりだった。彼の目標は、このパターンを最速で計算し、宇宙の真実に迫ることだった。

アレックスは、最先端の技術を駆使してこの問題に取り組んでいた。彼は最新のコンピュータを用い、AVX-512命令という強力なSIMD命令セットを活用して、計算速度を極限まで高めることを決意した。これにより、計算の効率性と精度を劇的に向上させることができると信じていた。

ある夜、アレックスはコンピュータに最終的なコードを投入し、計算を開始した。彼の心臓はドキドキと高鳴り、コンピュータのスクリーンに映るデータが急速に変化していくのを見守っていた。秒針が進むごとに、計算結果が鮮明な画像となって浮かび上がり、まるで宇宙の中の神秘的な銀河が描き出されているかのようだった。

「これが、最速の計算によって明らかになる宇宙の真実だ！」アレックスは興奮を抑えきれず、叫んだ。

計算が完了するやいなや、アレックスはその結果を眺めて驚愕した。計算時間は予想をはるかに下回る、信じられない速さで完了していた。AVX-512命令の力によって、彼のコンピュータは前例のない速度でデータを処理し、宇宙の神秘的なパターンを鮮やかに描き出していたのだ。

その瞬間、アレックスはただの科学者から、計算の魔術師へと変貌を遂げた。彼の業績は科学界を震撼させ、世界中の研究者たちが彼の成果に注目した。アレックスの名は「超速の伝説」として広まり、彼の計算結果は宇宙探査の新たな標準となった。

アレックスはその後も、数式の中に潜む美しさを探求し続けることを決意した。彼の挑戦はまだ始まったばかりで、彼の目の前には無限の可能性が広がっていることを知っていた。そして、彼の心の中には、AVX-512命令によって開かれた新たな扉への興奮と期待が常に燃えていた。

このショートストーリーが、AVX-512命令を駆使した計算のスピードとその影響を感じさせるものとなることを願っています。

AVX-512命令を有効にした最適化ありの実行: 計算時間: 0.20秒

CUDA GPU : 2.45秒

Google corab で実験しました。

AVX-512は、Advanced Vector Extensions (AVX) の512ビット版で、特にデータ並列処理において高い性能を提供します。これを利用することで、さらに高速化が期待できますが、AVX-512がサポートされているCPUが必要です。

AVX-512を有効にしてコンパイルする手順

コンパイラオプションの設定

-mavx512f オプションを使用してAVX-512命令を有効にします。また、-march=native を追加することで、使用しているCPUに最適なオプションが自動的に選ばれるようになります。

ポイント
AVX-512命令の使用: -mavx512f オプションを追加して、AVX-512命令を使用します。これにより、AVX-512命令セットを持つCPUでさらに高い性能が期待できます。

AVX-512命令を使用することで、特に数値計算やデータ並列処理が重いタスクに対して大きな性能向上が期待できますが、対応するハードウェアが必要です。

AVX-512命令を有効にした最適化ありの実行: 計算時間: 0.20秒

ChatGPT
AVX-512命令を有効にした最適化ありの実行で、計算時間が0.20秒となったのは非常に良い結果です。
これは、AVX-512命令がデータ並列処理に対して非常に効率的であることを示しています。

まとめ
AVX-512の効果: AVX-512命令を使用することで、計算が大幅に高速化されました。これは、AVX-512が提供する512ビットのベクトル処理能力が、計算負荷の高いタスクに対して特に有効であるためです。
他の最適化手法との比較: SIMD命令を使用しない最適化ありのC++コード（2.99秒）やCUDA GPU（2.45秒）と比べても、AVX-512を使用した最適化が優れた性能を発揮しています。

以下のコードは、AVX-512命令を有効にしてC++コードをコンパイルし、実行するためのものです：

import subprocess
import time
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

# C++のソースコード
cpp_code = """
#include <iostream>
#include <fstream>
#include <vector>
#include <complex>
#include <chrono>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace std::chrono;

int mandelbrot(const complex<double>& c, int max_iter) {
    complex<double> z = c;
    for (int n = 0; n < max_iter; ++n) {
        if (abs(z) > 2.0) {
            return n;
        }
        z = z * z + c;
    }
    return max_iter;
}

void plot_mandelbrot(double xmin, double xmax, double ymin, double ymax, int width, int height, int max_iter, const string& filename) {
    vector<uint8_t> img(width * height);
    double x_scale = (xmax - xmin) / (width - 1);
    double y_scale = (ymax - ymin) / (height - 1);

    for (int y = 0; y < height; ++y) {
        for (int x = 0; x < width; ++x) {
            double cr = xmin + x * x_scale;
            double ci = ymin + y * y_scale;
            complex<double> c(cr, ci);
            int color = mandelbrot(c, max_iter);
            img[y * width + x] = static_cast<uint8_t>(color % 256);
        }
    }

    cv::Mat image(height, width, CV_8UC1, img.data());
    cv::applyColorMap(image, image, cv::COLORMAP_INFERNO);
    cv::imwrite(filename, image);
}

int main() {
    auto start = high_resolution_clock::now();

    int width = 800, height = 600, max_iter = 256;
    vector<tuple<double, double, double, double>> positions = {
        {-2.0, 1.0, -1.5, 1.5},
        {-1.0, 0.5, -0.5, 0.5},
        {-1.5, -1.0, -1.0, 1.0},
        {-0.5, 0.5, -0.5, 0.5},
        {-2.5, -1.5, -1.5, 1.5}
    };

    for (size_t i = 0; i < positions.size(); ++i) {
        double xmin, xmax, ymin, ymax;
        tie(xmin, xmax, ymin, ymax) = positions[i];
        string filename = "mandelbrot_plot_" + to_string(i + 1) + ".png";
        plot_mandelbrot(xmin, xmax, ymin, ymax, width, height, max_iter, filename);
    }

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<milliseconds>(stop - start);
    cout << "計算時間: " << duration.count() << "ミリ秒" << endl;

    return 0;
}
"""

# C++ソースコードをファイルに書き込む
with open('mandelbrot.cpp', 'w') as file:
    file.write(cpp_code)

# OpenCVのパスを取得
pkg_config_output = subprocess.check_output(['pkg-config', '--cflags', '--libs', 'opencv4']).decode().strip()

# AVX-512命令を有効にして最適化ありでコンパイル
compile_command_opt_simd = f'g++ -O3 -march=native -mavx512f -o mandelbrot_opt_simd mandelbrot.cpp {pkg_config_output}'
try:
    subprocess.run(compile_command_opt_simd, shell=True, check=True)
except subprocess.CalledProcessError as e:
    print("コンパイルエラー:")
    print(e.output)

# AVX-512命令を有効にして最適化ありで実行し、処理時間を測定
print("\nAVX-512命令を有効にした最適化ありの実行:")
start_time = time.time()
result_opt_simd = subprocess.run(['./mandelbrot_opt_simd'], capture_output=True, text=True)
end_time = time.time()
print(result_opt_simd.stdout)
print(f'計算時間: {end_time - start_time:.2f}秒')

# 生成された画像を表示する
for i in range(1, 6):
    img = mpimg.imread(f'mandelbrot_plot_{i}.png')
    plt.figure(figsize=(10, 7))
    plt.imshow(img, cmap='inferno')
    plt.title(f'Mandelbrot Set - Plot {i}')
    plt.axis('off')
    plt.show()

参考。　Google corab環境。

no_optでの実行:
計算時間: 21456ミリ秒

計算時間: 21.93秒

optでの実行:
計算時間: 3266ミリ秒

計算時間: 3.38秒

avx2での実行:
計算時間: 5322ミリ秒

計算時間: 5.44秒

avx512での実行:

計算時間: 0.12秒

import subprocess
import time
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

# C++のソースコード
cpp_code = """
#include <iostream>
#include <fstream>
#include <vector>
#include <complex>
#include <chrono>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace std::chrono;

int mandelbrot(const complex<double>& c, int max_iter) {
    complex<double> z = c;
    for (int n = 0; n < max_iter; ++n) {
        if (abs(z) > 2.0) {
            return n;
        }
        z = z * z + c;
    }
    return max_iter;
}

void plot_mandelbrot(double xmin, double xmax, double ymin, double ymax, int width, int height, int max_iter, const string& filename) {
    vector<uint8_t> img(width * height);
    double x_scale = (xmax - xmin) / (width - 1);
    double y_scale = (ymax - ymin) / (height - 1);

    for (int y = 0; y < height; ++y) {
        for (int x = 0; x < width; ++x) {
            double cr = xmin + x * x_scale;
            double ci = ymin + y * y_scale;
            complex<double> c(cr, ci);
            int color = mandelbrot(c, max_iter);
            img[y * width + x] = static_cast<uint8_t>(color % 256);
        }
    }

    cv::Mat image(height, width, CV_8UC1, img.data());
    cv::applyColorMap(image, image, cv::COLORMAP_INFERNO);
    cv::imwrite(filename, image);
}

int main() {
    auto start = high_resolution_clock::now();

    int width = 800, height = 600, max_iter = 256;
    vector<tuple<double, double, double, double>> positions = {
        {-2.0, 1.0, -1.5, 1.5},
        {-1.0, 0.5, -0.5, 0.5},
        {-1.5, -1.0, -1.0, 1.0},
        {-0.5, 0.5, -0.5, 0.5},
        {-2.5, -1.5, -1.5, 1.5}
    };

    for (size_t i = 0; i < positions.size(); ++i) {
        double xmin, xmax, ymin, ymax;
        tie(xmin, xmax, ymin, ymax) = positions[i];
        string filename = "mandelbrot_plot_" + to_string(i + 1) + ".png";
        plot_mandelbrot(xmin, xmax, ymin, ymax, width, height, max_iter, filename);
    }

    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<milliseconds>(stop - start);
    cout << "計算時間: " << duration.count() << "ミリ秒" << endl;

    return 0;
}
"""

# C++ソースコードをファイルに書き込む
with open('mandelbrot.cpp', 'w') as file:
    file.write(cpp_code)

# OpenCVのパスを取得
pkg_config_output = subprocess.check_output(['pkg-config', '--cflags', '--libs', 'opencv4']).decode().strip()

# コンパイルオプションの定義
compile_options = [
    ('no_opt', 'g++ -o mandelbrot_no_opt mandelbrot.cpp {0}'),
    ('opt', 'g++ -O3 -o mandelbrot_opt mandelbrot.cpp {0}'),
    ('avx2', 'g++ -O3 -mavx2 -o mandelbrot_avx2 mandelbrot.cpp {0}'),
    ('avx512', 'g++ -O3 -mavx512f -o mandelbrot_avx512 mandelbrot.cpp {0}')
]

# コンパイルと実行時間の測定
execution_times = {}
for label, compile_cmd in compile_options:
    # コンパイル
    compile_command = compile_cmd.format(pkg_config_output)
    try:
        subprocess.run(compile_command, shell=True, check=True)
    except subprocess.CalledProcessError as e:
        print(f"コンパイルエラー({label}):")
        print(e.output)
        continue
    
    # 実行と時間測定
    print(f"\n{label}での実行:")
    start_time = time.time()
    result = subprocess.run([f'./mandelbrot_{label}'], capture_output=True, text=True)
    end_time = time.time()
    execution_times[label] = end_time - start_time
    print(result.stdout)
    print(f'計算時間: {execution_times[label]:.2f}秒')

# 実行時間のプロット
labels = list(execution_times.keys())
times = list(execution_times.values())

plt.bar(labels, times, color=['blue', 'green', 'red', 'purple'])
plt.xlabel('Optimization Type')
plt.ylabel('Execution Time (s)')
plt.title('Mandelbrot Set Calculation Time Comparison')
plt.show()

# 生成された画像を表示する
for i in range(1, 6):
    img = mpimg.imread(f'mandelbrot_plot_{i}.png')
    plt.figure(figsize=(10, 7))
    plt.imshow(img, cmap='inferno')
    plt.title(f'Mandelbrot Set - Plot {i}')
    plt.axis('off')
    plt.show()

CUDA GPU 処理時間: 2.32秒　　Google corab環境。(上記と同様の計算量で比較。)

# インストールコマンド
!pip install matplotlib numpy pycuda

import numpy as np
import matplotlib.pyplot as plt
import time
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule

# CUDAカーネルの定義
cuda_code = """
__global__ void mandelbrot(float *xmin, float *xmax, float *ymin, float *ymax, int *width, int *height, int *max_iter, int *img) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (x < *width && y < *height) {
        float r = *xmin + x * (*xmax - *xmin) / (*width - 1);
        float i = *ymin + y * (*ymax - *ymin) / (*height - 1);
        float cr = r;
        float ci = i;
        int iteration = 0;
        float zr = 0;
        float zi = 0;
        
        while (zr * zr + zi * zi <= 4.0f && iteration < *max_iter) {
            float temp = zr * zr - zi * zi + cr;
            zi = 2.0f * zr * zi + ci;
            zr = temp;
            iteration++;
        }
        
        img[y * (*width) + x] = iteration;
    }
}
"""

def plot_mandelbrot(xmin, xmax, ymin, ymax, width, height, max_iter):
    mod = SourceModule(cuda_code)
    mandelbrot = mod.get_function("mandelbrot")
    
    img = np.zeros((height, width), dtype=np.int32)
    
    # デバイスメモリの確保
    d_img = cuda.mem_alloc(img.nbytes)
    d_xmin = cuda.mem_alloc(np.float32(1).nbytes)
    d_xmax = cuda.mem_alloc(np.float32(1).nbytes)
    d_ymin = cuda.mem_alloc(np.float32(1).nbytes)
    d_ymax = cuda.mem_alloc(np.float32(1).nbytes)
    d_width = cuda.mem_alloc(np.int32(1).nbytes)
    d_height = cuda.mem_alloc(np.int32(1).nbytes)
    d_max_iter = cuda.mem_alloc(np.int32(1).nbytes)

    # ホストからデバイスへのデータ転送
    cuda.memcpy_htod(d_xmin, np.array([xmin], dtype=np.float32))
    cuda.memcpy_htod(d_xmax, np.array([xmax], dtype=np.float32))
    cuda.memcpy_htod(d_ymin, np.array([ymin], dtype=np.float32))
    cuda.memcpy_htod(d_ymax, np.array([ymax], dtype=np.float32))
    cuda.memcpy_htod(d_width, np.array([width], dtype=np.int32))
    cuda.memcpy_htod(d_height, np.array([height], dtype=np.int32))
    cuda.memcpy_htod(d_max_iter, np.array([max_iter], dtype=np.int32))

    # CUDAカーネルの実行
    block_size = (16, 16, 1)
    grid_size = (width // block_size[0] + 1, height // block_size[1] + 1)
    mandelbrot(d_xmin, d_xmax, d_ymin, d_ymax, d_width, d_height, d_max_iter, d_img, block=block_size, grid=grid_size)

    # デバイスからホストへのデータ転送
    cuda.memcpy_dtoh(img, d_img)
    
    return img

def main():
    # パラメータ設定
    width, height = 800, 600
    max_iter = 256
    positions = [
        (-2.0, 1.0, -1.5, 1.5),
        (-1.0, 0.5, -0.5, 0.5),
        (-1.5, -1.0, -1.0, 1.0),
        (-0.5, 0.5, -0.5, 0.5),
        (-2.5, -1.5, -1.5, 1.5)
    ]
    
    start_time = time.time()
    
    for idx, (xmin, xmax, ymin, ymax) in enumerate(positions):
        plt.figure(figsize=(10, 7))
        img = plot_mandelbrot(xmin, xmax, ymin, ymax, width, height, max_iter)
        plt.imshow(img, extent=(xmin, xmax, ymin, ymax), cmap='inferno')
        plt.colorbar()
        plt.title(f'Mandelbrot Set - Plot {idx+1}')
        plt.savefig(f'mandelbrot_plot_{idx+1}.png')
        plt.close()
    
    end_time = time.time()
    print(f'処理時間: {end_time - start_time:.2f}秒')

if __name__ == "__main__":
    main()