#include <complex>
#include <iostream>
#include <iterator>
#include <tuple>
#include <vector>

#include "../util/util.hpp"

// マンデルブロ集合の漸化式が何回で無限大に発散するか。発散しない場合は0を返します。
inline auto divergence_count(const std::complex<float>& c) noexcept {
    auto z = std::complex(0.0f, 0.0f);

    for (auto i = 0; i < 100; ++i) {
        if (std::isinf(z.real())) {
            return i;
        }

        z = std::pow(z, 2) + c;
    }

    return 0;
}

// real_minからreal_max、imag_minからimag_maxの範囲の複素数平面をreal_size×imag_sizeに区切って、それぞれがマンデルブロ集合に含まれるか確認します。
inline auto mandelbrot_set(float real_min, float real_max, std::size_t real_size, float imag_min, float imag_max, std::size_t imag_size) noexcept {
    auto result = std::vector<int>(); result.reserve(real_size * imag_size);

    for (const auto& imag: util::linspace<float>(imag_min, imag_max, imag_size)) {
        for (const auto& real: util::linspace<float>(real_min, real_max, real_size)) {
            result.emplace_back(divergence_count(std::complex(real, imag)));
        }
    }

    return std::make_tuple(result, real_size, imag_size);
}

// 適切なreal_sizeとimag_sizeを作成して、mandelbrot_set()を呼び出します。
inline auto mandelbrot_set(float real_min, float real_max, float imag_min, float imag_max) noexcept {
    auto [real_size, imag_size] = [&]() {  // ラムダ式を使うと変数のスコープを短くできて便利。構造化束縛を使うと複数の値をリターンできて便利。
        const auto real_diff = real_max - real_min;
        const auto imag_diff = imag_max - imag_min;

        return std::make_tuple(
            static_cast<std::size_t>(1000 * std::min(real_diff / imag_diff, 1.0f)),
            static_cast<std::size_t>(1000 * std::min(imag_diff / real_diff, 1.0f))
        );
    }();

    return mandelbrot_set(real_min, real_max, real_size, imag_min, imag_max, imag_size);で、
}

// メイン・ルーチン。
int main(int argc, char** argv) {
    // 時間計速。util::timeitは自作の処理時間計測ルーチン。./util/util.hppを参照してください。
    util::timeit([&]() {
        mandelbrot_set(-2.0f, 2.0f, -2.0f, 2.0f);
    });

    // プログラムが正しいかを確認するため、マンデルブロ集合を出力します。
    const auto [v, w, h] = mandelbrot_set(-2.0f, 2.0f, -2.0f, 2.0f);

    {
        auto it = std::begin(v);

        for (auto i = 0_z; i < h; ++i, it += w) {
            std::copy(it, it + w, std::ostream_iterator<float>(std::cout, "\t"));
            std::cout << std::endl;
        }
    }

    return 0;
}

#include <complex>
#include <iostream>
#include <iterator>
#include <tuple>
#include <vector>

#include <cuda.h>
#include <thrust/complex.h>

#include "../util/cudautil.hpp"
#include "../util/util.hpp"

__device__  // GPU上で実行される関数には__device__を付けてください。
inline auto divergence_count(const thrust::complex<float>& c) noexcept {
    auto z = thrust::complex(0.0f, 0.0f);  // __device__関数ではライブラリの呼び出しに制限があって、std::complexが使えません。Thrustのcomplexで代用しました。

    for (auto i = 0; i < 100; ++i) {
        if (std::isinf(z.real())) {
            return i;
        }

        z = z * z + c;  // thrust::complexではstd::pow()を使用できなかったので、z * zに書き換えました。
    }

    return 0;
}

__device__
inline auto linspace(float start, float stop, std::size_t size) noexcept {  // util::linspace()は__device__ではないので、__device__バージョンを作成しました。
    auto result = new float[size];  // std::vectorは使用できなかったので、生ポインターで。

    const auto delta = (stop - start) / static_cast<float>(size - 1);

    for (auto i = static_cast<std::size_t>(0); i < size; ++i) {
        result[i] = start + delta * static_cast<float>(i);
    }

    return result;
}

__global__  // GPU上で実行され、かつ、CPUから呼び出される関数には__global__をつけてください。
void mandelbrot_set(float real_min, float real_max, std::size_t real_size, float imag_min, float imag_max, std::size_t imag_size, int* result) {  // __global__の戻り値はvoidです。また、noexceptも付加できません。
    const auto imags = linspace(imag_min, imag_max, imag_size);  // あとでdeleteしなければならないので、変数に保存しておきます。
    const auto reals = linspace(real_min, real_max, real_size);

    for (auto i = 0; i < imag_size; ++i) {  // ポインターだと拡張forが使えないので、普通のforを使います。
        for (auto j = 0; j < real_size; ++j) {
            result[i * real_size + j] = divergence_count(thrust::complex(reals[j], imags[i]));
        }
    }

    delete reals;  // linspaceはnewしたポインターを返すので、自前でdeleteします。
    delete imags;
}

inline auto mandelbrot_set(float real_min, float real_max, float imag_min, float imag_max) noexcept {
    const auto [real_size, imag_size] = [&]() {
        const auto real_diff = real_max - real_min;
        const auto imag_diff = imag_max - imag_min;

        return std::make_tuple(
            static_cast<std::size_t>(1000 * std::min(real_diff / imag_diff, 1.0f)),
            static_cast<std::size_t>(1000 * std::min(imag_diff / real_diff, 1.0f))
        );
    }();

    auto result = std::vector<int>(real_size * imag_size);

    int* result_device;  // GPUのメモリとCPUのメモリは異なるので、GPU用のメモリを用意します。
    cuda_check(cudaMalloc(&result_device, sizeof(int) * real_size * imag_size));  // GPUのメモリの確保はcudaMallocで。エラー・チェックは戻り値でやるので、cudautil::cuda_check()という関数を自作しました。

    mandelbrot_set<<<1, 1>>>(real_min, real_max, real_size, imag_min, imag_max, imag_size, result_device);  // __global__の呼び出しには、<<<block, grid>>>で並列度を指定します。今回は1並列（並列実行しない）でやります。
    cuda_check(cudaGetLastError());  // __global__の呼び出しでエラーが発生してないかのチェック。

    cuda_check(cudaMemcpy(result.data(), result_device, sizeof(int) * real_size * imag_size, cudaMemcpyDeviceToHost));  // GPUのメモリーからCPUのメモリーに結果をコピーします。

    cuda_check(cudaFree(result_device));  // GPUのメモリ開放はcudaFreeでやります。

    return std::make_tuple(result, real_size, imag_size);
}

int main(int argc, char** argv) {
    first_cudaMalloc_is_too_slow();  // 一回目のcudaMallocは遅いので、タイム計測に影響しないように一回cudaMallocしておきます。
    cuda_check(cudaDeviceSynchronize());  // CUDAの関数呼び出しは非同期なので、タイム計測に影響しないようにGPU上の処理が終わるまで待ちます。

    util::timeit([&]() {
        mandelbrot_set(-2.0f, 2.0f, -2.0f, 2.0f);
    });

    const auto [v, w, h] = mandelbrot_set(-2.0f, 2.0f, -2.0f, 2.0f);

    {
        auto it = std::begin(v);

        for (auto i = 0_z; i < h; ++i, it += w) {
            std::copy(it, it + w, std::ostream_iterator<float>(std::cout, "\t"));
            std::cout << std::endl;
        }
    }

    return 0;
}

#include <cstdlib>

#include <cuda.h>

inline auto cuda_check(const cudaError_t& error) noexcept {
    if (error == cudaSuccess) {
        return;
    }

    std::cerr << cudaGetErrorString(error) << std::endl;
    std::exit(1);
}

#include <iostream>
#include <iterator>
#include <tuple>
#include <vector>

#include <cuda.h>
#include <thrust/complex.h>

#include "../util/cudautil.hpp"
#include "../util/util.hpp"

// mandelbrot-cuda-ngと同じ。標準ライブラリを使用しない形に書き換えただけです。
__device__
inline auto divergence_count(const thrust::complex<float>& c) noexcept {
    auto z = thrust::complex(0.0f, 0.0f);

    for (auto i = 0; i < 100; ++i) {
        if (std::isinf(z.real())) {
            return i;
        }

        z = z * z + c;
    }

    return 0;
}

// この関数をGPU上で並列で実行します。
__global__
void mandelbrot_set(float real_min, float real_max, std::size_t real_size, float imag_min, float imag_max, std::size_t imag_size, int* result) {
    const auto i = blockIdx.y * blockDim.y + threadIdx.y;  // 処理対象を取得します。詳細は本文を参照してください。
    const auto j = blockIdx.x * blockDim.x + threadIdx.x;  // 同上。

    if (i >= imag_size || j >= real_size) {  // 範囲チェック。iがimag_sizeを超えることはないのですけど、パターン化しておいたほうがバグが減るので。
        return;
    }

    const auto imag = imag_min + (imag_max - imag_min) / static_cast<float>(imag_size - 1) * static_cast<float>(i);  // linspace()の処理を移植。スレッドのIDから、対象となる複素数を計算します。
    const auto real = real_min + (real_max - real_min) / static_cast<float>(real_size - 1) * static_cast<float>(j);  // 同上。

    result[i * real_size + j] = divergence_count(thrust::complex(real, imag));
}

// GPU上で実行する関数を、この関数から呼び出します。
inline auto mandelbrot_set(float real_min, float real_max, std::size_t real_size, float imag_min, float imag_max, std::size_t imag_size) noexcept {
    auto result = std::vector<int>(real_size * imag_size);

    int* result_device;
    cuda_check(cudaMalloc(&result_device, sizeof(int) * real_size * imag_size));

    // 並列でmandelbrot_setを実行します。gridとblockについては本文を参照してください。
    const auto grid  = dim3((real_size + 256 - 1) / 256, imag_size);
    const auto block = dim3(256, 1);
    mandelbrot_set<<<grid, block>>>(real_min, real_max, real_size, imag_min, imag_max, imag_size, result_device);
    cuda_check(cudaGetLastError());

    cuda_check(cudaMemcpy(result.data(), result_device, sizeof(int) * real_size * imag_size, cudaMemcpyDeviceToHost));

    cuda_check(cudaFree(result_device));

    return std::make_tuple(result, real_size, imag_size);
}

// mandelbrotと同じ。
inline auto mandelbrot_set(float real_min, float real_max, float imag_min, float imag_max) noexcept {
    const auto [real_size, imag_size] = [&]() {
        const auto real_diff = real_max - real_min;
        const auto imag_diff = imag_max - imag_min;

        return std::make_tuple(
            static_cast<std::size_t>(1000 * std::min(real_diff / imag_diff, 1.0f)),
            static_cast<std::size_t>(1000 * std::min(imag_diff / real_diff, 1.0f))
        );
    }();

    return mandelbrot_set(real_min, real_max, real_size, imag_min, imag_max, imag_size);
}

// mandelbrot-cuda-ngと同じ。
int main(int argc, char** argv) {
    first_cudaMalloc_is_too_slow();
    cuda_check(cudaDeviceSynchronize());

    util::timeit([&]() {
        mandelbrot_set(-2.0f, 2.0f, -2.0f, 2.0f);
    });

    const auto [v, w, h] = mandelbrot_set(-2.0f, 2.0f, -2.0f, 2.0f);

    {
        auto it = std::begin(v);

        for (auto i = 0_z; i < h; ++i, it += w) {
            std::copy(it, it + w, std::ostream_iterator<float>(std::cout, "\t"));
            std::cout << std::endl;
        }
    }

    return 0;
}

const auto grid  = dim3((real_size + 256 - 1) / 256, imag_size);
const auto block = dim3(256, 1);
mandelbrot_set<<<grid, block>>>(real_min, real_max, real_size, imag_min, imag_max, imag_size, result_device);

const auto i = blockIdx.y * blockDim.y + threadIdx.y;
const auto j = blockIdx.x * blockDim.x + threadIdx.x;

#include <iostream>
#include <iterator>
#include <tuple>
#include <vector>

#include <cuda.h>
#include <thrust/complex.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/transform.h>

#include "../util/cudautil.hpp"
#include "../util/util.hpp"

// mandelbrot-cudaと同じ。標準ライブラリを使用しない形に書き換えただけです。
__device__
inline auto divergence_count(const thrust::complex<float>& c) noexcept {
    auto z = thrust::complex(0.0f, 0.0f);

    for (auto i = 0; i < 100; ++i) {
        if (std::isinf(z.real())) {
            return i;
        }

        z = z * z + c;
    }

    return 0;
}

// この関数の内部で、Thrustを使用して並列処理します。Thrustのおかげで__global__は不要になりました。
inline std::tuple<thrust::host_vector<int>, std::size_t, std::size_t> mandelbrot_set(float real_min, float real_max, std::size_t real_size, float imag_min, float imag_max, std::size_t imag_size) noexcept {  // __device__なラムダ式を含む関数では、戻り値にautoを使用できません……。
    auto result_device = thrust::device_vector<int>(real_size * imag_size);  // thrust::device_vectorは、GPUメモリを使うvector。
    auto it = thrust::counting_iterator<std::size_t>(0);  // thrust::counting_iteratorは数字をカウント・アップしてくれるイテレーター。最初の要素は0、次の要素は1になります。
    thrust::transform(it, it + real_size * imag_size, std::begin(result_device),  // real_size * imag_size個の要素をtransformして、result_deviceに代入します。
        [=] __device__ (const auto& i) {  // __device__なラムダ式は、[&]は使用できません。[=]を使用してください。
            const auto imag = imag_min + (imag_max - imag_min) / static_cast<float>(imag_size - 1) * static_cast<float>(i / imag_size);  // インデックスから複素数の虚数部を計算します。
            const auto real = real_min + (real_max - real_min) / static_cast<float>(real_size - 1) * static_cast<float>(i % imag_size);  // インデックスから複素数の実数部を計算します。

            return divergence_count(thrust::complex(real, imag));
        }
    );

    return std::make_tuple(thrust::host_vector<int>(result_device), real_size, imag_size);  // thrust::host_vectorにコピーしてリターンします。
}

// mandelbrot-cudaと同じ。
inline auto mandelbrot_set(float real_min, float real_max, float imag_min, float imag_max) {
    const auto [real_size, imag_size] = [&]() {
        const auto real_diff = real_max - real_min;
        const auto imag_diff = imag_max - imag_min;

        return std::make_tuple(
            static_cast<std::size_t>(1000 * std::min(real_diff / imag_diff, 1.0f)),
            static_cast<std::size_t>(1000 * std::min(imag_diff / real_diff, 1.0f))
        );
    }();

    return mandelbrot_set(real_min, real_max, real_size, imag_min, imag_max, imag_size);
}

// mandelbrot-cudaと同じ。
int main(int argc, char** argv) {
    first_cudaMalloc_is_too_slow();
    cuda_check(cudaDeviceSynchronize());

    util::timeit([&]() {
        mandelbrot_set(-2.0f, 2.0f, -2.0f, 2.0f);

        cuda_check(cudaDeviceSynchronize());
    });

    const auto [v, w, h] = mandelbrot_set(-2.0f, 2.0f, -2.0f, 2.0f);

    {
        auto it = std::begin(v);

        for (auto i = 0_z; i < h; ++i, it += w) {
            std::copy(it, it + w, std::ostream_iterator<float>(std::cout, "\t"));
            std::cout << std::endl;
        }
    }

    return 0;
}

#include <iostream>
#include <random>

#include "../util/util.hpp"

// 円の内側にあるかを判断します。
inline auto is_in_circle(float x, float y) noexcept {
    return std::pow(x, 2) + std::pow(y, 2) <= 1.0f;  // 原点は(0, 0)。乱数を0〜1にしたので半径は1となるので、平方根を求める処理は不要。
}

// モンテカルロ法で円周率を求めます。
inline auto monte_carlo_pi(int n, unsigned int seed) noexcept {
    auto rng = std::mt19937(seed);  // 乱数ジェネレーター。メルセンヌ・ツイスタを使用します。
    auto distribution = std::uniform_real_distribution(0.0f, 1.0f);  // 乱数は、0〜1の一様分布にします。マイナスの値がないので、円の右上の1/4の範囲を調べることになります。

    auto c = 0;  // 円の内側と判定された数。c / nは、円の面積のpi * r * rの1/4と等しくなるはす。

    for (auto i = 0; i < n; ++i) {  // n回繰り返して……
        if (is_in_circle(distribution(rng), distribution(rng))) {  // ランダムな点が円の内側にあるなら……
            c++;  // cをインクリメントします。
        }
    }

    return 4.0f * static_cast<float>(c) / static_cast<float>(n);  // このコードでの半径は1なので、c / n = pi * 1 * 1 / 4。だから4 * c / n = piになるはず。
}

// メイン・ルーチン
int main(int argc, char** argv) {
    util::timeit([]() {  // 処理時間を計測します。
        monte_carlo_pi(100'000'000, 0u);
    });

    std::cout << monte_carlo_pi(100'000'000, 0ul) << std::endl;  // 正しいかを確認するために、処理結果を出力します。

    return 0;
}

#include <iostream>
#include <random>

#include <cuda.h>
#include <curand.h>
#include <thrust/count.h>
#include <thrust/device_vector.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/tuple.h>

#include "../util/cudautil.hpp"
#include "../util/util.hpp"

// __device__を追加したただけ。他はmonte-carloと同じ。
__device__
inline auto is_in_circle(float x, float y) noexcept {
    return std::pow(x, 2) + std::pow(y, 2) <= 1.0f;
}

inline float monte_carlo_pi(int n, unsigned long long seed) noexcept {
    auto rng = curandGenerator_t();  // 乱数ジェネレーター。
    curand_check(curandCreateGenerator(&rng, CURAND_RNG_PSEUDO_MTGP32));  // 乱数ジェネレーターを初期化します。MTGP32は、GPUに最適化されたメルセンヌ・ツイスタ。
    curand_check(curandSetPseudoRandomGeneratorSeed(rng, seed));  // シードを設定します。

    auto xs = thrust::device_vector<float>(n);  // メモリ管理で楽をしたいのでthrust::device_vectorを使用します。
    curand_check(curandGenerateUniform(rng, xs.data().get(), n));  // 一様乱数を生成します。thrust::device_vectorのポインターはdata().get()で取得できます。

    auto ys = thrust::device_vector<float>(n);  // yについても、同様に生成します。
    curand_check(curandGenerateUniform(rng, ys.data().get(), n));

    curand_check(curandDestroyGenerator(rng));  // 乱数ジェネレーターを破棄します。

    auto it = thrust::make_zip_iterator(thrust::make_tuple(std::begin(xs), std::begin(ys)));  // xとyで2つの乱数が必要なので、thrust::make_zip_iteratorで一つにまとめます。
    const auto c = thrust::count_if(it, it + n,  // thrust::count_ifで、ラムダ式が真となった場合の数を数えます。
        [=] __device__ (const auto& p) {
            return is_in_circle(thrust::get<0>(p), thrust::get<1>(p));  // make_zip_iteratorはthrust::tupleを要素に持つので、get<0>()とget<1>()でタプルの要素を取得します。
        }
    );

    return 4.0f * static_cast<float>(c) / static_cast<float>(n);
}

int main(int argc, char** argv) {
    first_cudaMalloc_is_too_slow();  // 一回目のcudaMallocは遅いので、タイム計測に影響しないように一回cudaMallocしておきます。
    cuda_check(cudaDeviceSynchronize());  // CUDAの関数呼び出しは非同期なので、タイム計測に影響しないようにGPU上の処理が終わるまで待ちます。

    util::timeit([]() {
        monte_carlo_pi(100'000'000, 0ull);
    });

    std::cout << monte_carlo_pi(100'000'000, 0ul) << std::endl;

    return 0;
}

#include <iostream>
#include <vector>

#include "../util/util.hpp"

const std::size_t M = 1001;
const std::size_t N = 1002;
const std::size_t K = 1003;

// size_m行size_k列の行列marix_aと、size_k行size_n列のmartix_bを掛け算して、matrix_cに結果を入れます。matrix_cは容量を確保した上で、0に初期化しておいてください。
inline auto matmul(const std::vector<float>& matrix_a, const std::vector<float>& matrix_b, std::vector<float>& matrix_c, std::size_t size_m, std::size_t size_n, std::size_t size_k) noexcept {
    for (auto i = 0_z; i < size_m; ++i) {
        for (auto j = 0_z; j < size_n; ++j) {
            // 行列の掛け算は、matrix_aの横に並ぶ要素とmatrix_bの縦に並ぶ要素を掛け算して、足し合わせるらしい。
            for (auto k = 0_z; k < size_k; ++k) {
                matrix_c[i * size_n + j] += matrix_a[i * size_k + k] * matrix_b[k * size_n + j];
            }
        }
    }
}

int main(int argc, char** argv) {
    const auto a = util::linspace(0.0f, 1.0f, M * K);  // M行K列の行列。
    const auto b = util::linspace(0.0f, 1.0f, K * N);  // K行M列の行列。

    // 処理時間を計測します。
    util::timeit([&]() {
        auto c = std::vector<float>(M * N, 0.0f);  // 掛け算の結果を入れるM行N列の行列。0で初期化しておきます。
        matmul(a, b, c, M, N, K);  // aとbを掛け算して、結果をcに格納します。
    });

    auto c = std::vector<float>(M * N, 0.0f);
    matmul(a, b, c, M, N, K);

    ...

    return 0;
}

for (auto k = 0_z; k < size_k; ++k) {
    matrix_c[i * size_n + j] += matrix_a[i * size_k + k] * matrix_b[k * size_n + j];
}

inline auto matmul_ng(const std::vector<float>& matrix_a, const std::vector<float>& matrix_b, std::vector<float>& matrix_c, std::size_t size_m, std::size_t size_n, std::size_t size_k) noexcept {
    for (auto i = 0_z; i < size_m; ++i) {
        for (auto k = 0_z; k < size_k; ++k) {  // この行と……
            for (auto j = 0_z; j < size_n; ++j) {  // この行を入れ替えた。
                matrix_c[i * size_n + j] += matrix_a[i * size_k + k] * matrix_b[k * size_n + j];  // 最も内側のループ実行時にアクセスするメモリ・アドレスが隣り合っているので、キャッシュ・メモリが使用されて効率が良くなる。
            }
        }
    }
}

#include <iostream>
#include <vector>

#include <cuda.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>

#include "../util/cudautil.hpp"
#include "../util/util.hpp"

const std::size_t M = 1001;
const std::size_t N = 1002;
const std::size_t K = 1003;

// 並列で実行される行列の掛け算。
__global__
void matmul(const float* matrix_a, const float* matrix_b, float* matrix_c, std::size_t size_m, std::size_t size_n, std::size_t size_k) {
    const auto i = blockIdx.y * blockDim.y + threadIdx.y;  // スレッドで処理すべき対象を取得します。
    const auto j = blockIdx.x * blockDim.x + threadIdx.x;  // スレッドで処理すべき対象を取得します。

    if (i >= size_m || j >= size_n) {  // 範囲チェック。
        return;
    }

    // matrix_aの横に並ぶ要素とmatrix_bの縦に並ぶ要素を掛け算して、足し合わせます。
    for (auto k = static_cast<std::size_t>(0); k < size_k; ++k) {
        matrix_c[i * size_n + j] += matrix_a[i * size_k + k] * matrix_b[k * size_n + j];
    }
}

int main(int argc, char** argv) {
    const auto a = thrust::host_vector<float>(util::linspace(0.0f, 1.0f, M * K));  // M行K列の行列。使用するメモリはCPUのメモリ。
    const auto a_device = thrust::device_vector<float>(a);  // M行K列の行列。使用するメモリはGPUのメモリ。

    const auto b = thrust::host_vector<float>(util::linspace(0.0f, 1.0f, K * N));  // K行N列の行列。使用するメモリはCPUのメモリ。
    const auto b_device = thrust::device_vector<float>(b);  // K行N列の行列。使用するメモリがGPUのメモリ。

    cuda_check(cudaDeviceSynchronize());  // CUDAの関数呼び出しは非同期なので、タイム計測に影響しないようにGPU上の処理が終わるまで待ちます。

    util::timeit([&]() {
        auto c_device = thrust::device_vector<float>(M * N, 0.0f);  // 掛け算の結果を入れるM行N列の行列。使用するメモリはGPUのメモリ。0で初期化しておきます。

        matmul<<<dim3((M + 16 - 1) / 16, (N + 16 - 1) / 16), dim3(16, 16)>>>(a_device.data().get(), b_device.data().get(), c_device.data().get(), M, N, K);  // a_deviceとb_deviceを掛け算して、結果をc_deviceに格納します。
        cuda_check(cudaGetLastError());

        cuda_check(cudaDeviceSynchronize());
    });

    auto c_device = thrust::device_vector<float>(M * N, 0.0f);

    matmul<<<dim3((M + 16 - 1) / 16, (N + 16 - 1) / 16), dim3(16, 16)>>>(a_device.data().get(), b_device.data().get(), c_device.data().get(), M, N, K);
    cuda_check(cudaGetLastError());

    const auto c = thrust::host_vector<float>(c_device);

    ...

    return 0;
}

__global__
void matmul(const float* matrix_a, const float* matrix_b, float* matrix_c, std::size_t size_m, std::size_t size_n, std::size_t size_k) {
    const auto i = blockIdx.y * blockDim.y + threadIdx.y;
    const auto k = blockIdx.x * blockDim.x + threadIdx.x;  // jからkに変更しました。

    if (i >= size_m || k >= size_k) {
        return;
    }

    for (auto j = static_cast<std::size_t>(0); j < size_n; ++j) {  // kのループをjのループに変更しました。
        matrix_c[i * size_n + j] += matrix_a[i * size_k + k] * matrix_b[k * size_n + j];
    }
}

matmul<<<dim3((M + 16 - 1) / 16, (K + 16 - 1) / 16), dim3(16, 16)>>>(a_device.data().get(), b_device.data().get(), c_device.data().get(), M, N, K);  // JではなくKで並列化するように変更しました。

#include <iostream>
#include <vector>

#include <cublas_v2.h>
#include <cuda.h>
#include <thrust/device_vector.h>

#include "../util/cudautil.hpp"
#include "../util/util.hpp"

const std::size_t M = 1001;
const std::size_t N = 1002;
const std::size_t K = 1003;

// 行列を転置します。M行N列の行列を、N行M列の行列に変換します。
template <typename T>
inline auto transpose(const std::vector<T>& matrix, std::size_t row_size, std::size_t col_size) noexcept {
    auto result = std::vector<T>(col_size * row_size);

    for (auto i = 0_z; i < col_size; ++i) {
        for (auto j = 0_z; j < row_size; ++j) {
            result[i * row_size + j] = matrix[j * col_size + i];  // 仮実装ということで、効率は無視させてください……。
        }
    }

    return result;
}

int main(int argc, char** argv) {
    cublasHandle_t handle;  // cuBLASのハンドル。
    cublas_check(cublasCreate(&handle));  // ハンドルを初期化します。

    const auto a = transpose(util::linspace(0.0f, 1.0f, M * K), M, K);  // cuBLASは列優先なので、今回はtransposeして対応しました。
    auto a_device = thrust::device_vector<float>(M * K);  // GPUのメモリを確保します。
    cublas_check(cublasSetMatrix(M, K, sizeof(float), a.data(), M, a_device.data().get(), M));  // cublasSetMatrixでデータを転送します。

    const auto b = transpose(util::linspace(0.0f, 1.0f, K * N), K, N);
    auto b_device = thrust::device_vector<float>(K * N);
    cublas_check(cublasSetMatrix(K, N, sizeof(float), b.data(), K, b_device.data().get(), K));

    cuda_check(cudaDeviceSynchronize());

    util::timeit([&]() {
        auto c_device = thrust::device_vector<float>(M * N);

        const auto alpha = 1.0f;  // cuBLASの引数はなんでかアドレス渡しなので、引数に必要な変数を作成します。
        const auto beta  = 0.0f;

        cublas_check(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, M, N, K, &alpha, a_device.data().get(), M, b_device.data().get(), K, &beta, c_device.data().get(), M));  // cublasSgemm()を呼び出します。

        cuda_check(cudaDeviceSynchronize());
    });

    auto c_device = thrust::device_vector<float>(M * N);

    const auto alpha = 1.0f;
    const auto beta  = 0.0f;

    cublas_check(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, M, N, K, &alpha, a_device.data().get(), M, b_device.data().get(), K, &beta, c_device.data().get(), M));

    auto c = std::vector<float>(N * M, 0.0f);
    cublas_check(cublasGetMatrix(M, N, sizeof(float), c_device.data().get(), M, c.data(), M));

    cublas_check(cublasDestroy(handle));  // ハンドルを破棄します。

    ...

    return 0;
}

// 普通のC++で行列を作る。
auto m = std::vector<float>(y_size * x_size);

for (auto y = 0; y < y_size; ++y) {  // まずは行のループ
    for (auto x = 0; x < x_size; ++x) {  // その中で列のループ
        m[y * x_size + x] = ...
    }
}

// 列優先の行列を作る。
auto m = std::vector<float>(x_size * y_size);

for (auto x = 0; x < x_size; ++x) {  // まずは列のループ
    for (auto y = 0; y < y_size; ++y) {  // その中で行のループ
        m[x * y_size + y] = ...
    }
}