Вероятно, вам уже попадались подобные руководства по CUDA: хрестоматийный пример «Hello World», в котором перемешан код для ЦП и графического процессора. Всё это сложено в один гетерогенный файл с исходниками на CUDA C++, а для запуска ядра применяется синтаксис NVCC с тройными угловыми скобками <<<>>>, который уже стал культовым:

#include <cuda_runtime.h>
#include <stdio.h>

__global__ void kernel() {
    printf("Hello World from block %d, thread %d\n", blockIdx.x, threadIdx.x);
}

int main() {
    kernel<<<1, 1>>>(); // Возвращает `void`?!     
    return cudaDeviceSynchronize() == cudaSuccess ? 0 : -1;
}

Время идёт, а такой паттерн по-прежнему попадается мне в продакшен-коде. Признаюсь, кое-где он всплывает и в моих любительских проектах — раз, два, три. Но это не лучшая идея, полагаться в серьёзном коде на запуск ядра через тройные угловые скобки. В таком случае программа не возвращает коды ошибок, поэтому может показаться обманчиво простой. Ниже вас ждут примерно 25 килобайт текста, в которых мы обсудим не самые корявые способы запуска ядер.

Основы и корректность

Вышеприведённый код скомпилируется, после чего выполнит ожидаемый вывод:

$ nvcc -o hello_world hello_world.cu && ./hello_world
> Hello World from block 0, thread 0

В каком-то смысле его уже можно считать «корректным».

Но в наше время обычны системы, в которых предусматривается не менее 8 GPU на плату HGX, поэтому, естественно, хочется реализовать некоторый параллелизм — чтобы гарантировать, что будут запускаться ядра на каждом из них.

Опуская базовые аспекты эксплуатации, отмечу, что каждый узел DGX H100 оснащается двумя могучими ядрами ЦП, которые без перерыва жонглируют сложными графами выполнения, одновременно передавая на вход и на выход сотни гигабайт. Всё это — асинхронно.

В такой системе очень многое может сбиться, так что давайте сформулируем несколько основополагающих правил оркестрации ядер GPU:

• Запуск каждого ядра сопряжён с серьёзной задержкой, поэтому такие операции должны выполняться асинхронно.

• Работу в пределах потоков следует раздавать явно.

• Вызовы API CUDA и пуски ядра должны сопровождаться надёжной проверкой на случай возможных ошибок.

Вот как интегрировать потоки CUDA — ничего сверхъестественного:

#include <cuda_runtime.h>
#include <stdio.h>

__global__ void kernel() {
    extern __shared__ char shared_buffer[];
    printf("Hello World from block %d, thread %d\n", blockIdx.x, threadIdx.x);
}

int main() {
    cudaStream_t stream;
    cudaStreamCreate(&stream);
    uint shared_memory_size = 0;
    kernel<<<1, 1, shared_memory_size, stream>>>(); // 4 аргумента, а не 2
    cudaStreamSynchronize(stream);
    cudaStreamDestroy(stream);
    return 0;
}

Обратите внимание: ядро запускается с четырьмя аргументами.

Вот более аккуратная версия с явной обработкой ошибок:

#include <cuda_runtime.h>
#include <stdio.h>

__global__ void kernel() {
    extern __shared__ char shared_buffer[];
    printf("Hello World from block %d, thread %d\n", blockIdx.x, threadIdx.x);
}

int main() {
    cudaStream_t stream;
    cudaError_t err = cudaStreamCreate(&stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to create stream: %s\n", cudaGetErrorString(err));
        return -1;
    }
    uint shared_memory_size = 1 << 30; // 1 ГБ – это много, но в демонстрационных целях такая величина подобрана специально
    kernel<<<1, 1, shared_memory_size, stream>>>();
    err = cudaStreamSynchronize(stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to synchronize stream: %s\n", cudaGetErrorString(err));
        cudaStreamDestroy(stream);
        return -1;
    }
    cudaStreamDestroy(stream);
    return 0;
}

Именно на этом, как правило, заканчиваются туториалы — и незаметно начинаются проблемы. Попробуйте запустить вот это:

$ nvcc -o hello_world hello_world.cu && ./hello_world

Никакого вывода, никакого сообщения об ошибке, вообще ничего. Мы не можем вытянуть ошибку из потока, поскольку отказ происходит при отправке данных, а не при выполнении.

API среды выполнения CUDA

Тройные угловые скобки NVCC — это синтаксический сахар, которым присыпан API среды выполнения CUDA, а этот интерфейс, в свою очередь, обёртывает более низкоуровневый API драйверов CUDA. Для эффективного отлова ошибок необходимо задействовать CUDA Driver Execution Control API:

#include <cuda_runtime.h>
#include <stdio.h>

__global__ void kernel() {
    extern __shared__ char shared_buffer[];
    printf("Hello World from block %d, thread %d\n", blockIdx.x, threadIdx.x);
}

int main() {
    cudaStream_t stream;
    cudaError_t err;
    err = cudaStreamCreate(&stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to create stream: %s\n", cudaGetErrorString(err));
        return -1;
    }

    dim3 grid(1);
    dim3 block(1);
    size_t shared_memory_size = 1 << 30; // 1 ГБ
    void *kernel_args[] = {};
    err = cudaLaunchKernel((void *)kernel, grid, block, kernel_args, shared_memory_size, stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to launch kernel: %s\n", cudaGetErrorString(err));
        cudaStreamDestroy(stream);
        return -1;
    }

    err = cudaStreamSynchronize(stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Kernel execution failed: %s\n", cudaGetErrorString(err));
        cudaStreamDestroy(stream);
        return -1;
    }
    err = cudaStreamDestroy(stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to destroy stream: %s\n", cudaGetErrorString(err));
        return -1;
    }
    return 0;
}

Скомпилируем и выполним:

$ nvcc -o hello_world hello_world.cu && ./hello_world
> Failed to launch kernel: invalid argument

Эта ошибка ожидаема, поскольку ядро невозможно даже отправить — из-за абсурдного запроса к памяти. Но проблема с этим API в том, что нам требуется найти иной способ, чтобы передавать аргументы ядру. Вот как мы передавали бы ядру массивы и скаляры:

#include <cuda_runtime.h>
#include <stdio.h>

__global__ void kernel(float *amount, size_t count, int power) {
    size_t idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx > count) return;
    amount[idx] = amount[idx] * scalbln(1.0, power); // Пример встроенной функции CUDA ;)
}

int main() {
    cudaError_t err;
    size_t num_elements = 1024;
    int integral_power = -2;
    double *data;

    // Выделяем объединённую память
    err = cudaMallocManaged(&data, num_elements * sizeof(double));
    if (err != cudaSuccess) {
        fprintf(stderr, "cudaMallocManaged failed: %s\n", cudaGetErrorString(err));
        return -1;
    }

    // Инициализируем данные
    for (size_t i = 0; i < num_elements; ++i) data[i] = (double)i;

    // Создаём поток CUDA 
    cudaStream_t stream;
    err = cudaStreamCreate(&stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to create stream: %s\n", cudaGetErrorString(err));
        cudaFree(data);
        return -1;
    }

    // Определяем параметры пуска ядра
    dim3 grid((num_elements + 255) / 256);
    dim3 block(256);
    void *kernel_args[] = {
        (void *)&data,
        (void *)&num_elements,
        (void *)&integral_power,
    };

    // Запускаем ядро
    err = cudaLaunchKernel((void *)kernel, grid, block, kernel_args, 0, stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to launch kernel: %s\n", cudaGetErrorString(err));
        cudaStreamDestroy(stream);
        cudaFree(data);
        return -1;
    }

    // Синхронизируем поток
    err = cudaStreamSynchronize(stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Kernel execution failed: %s\n", cudaGetErrorString(err));
        cudaStreamDestroy(stream);
        cudaFree(data);
        return -1;
    }

    // Выводим результаты на экран
    for (size_t i = 0; i < 5; ++i) printf("data[%zu] = %f\n", i, data[i]);
    cudaStreamDestroy(stream);
    cudaFree(data);
    return 0;
}

Я воспользовался объединённой памятью, чтобы упростить пример. Мы не обязаны явно выделять 2 буфера одновременно в ЦП и графическом процессоре и копировать данные между ними. Драйвер держит копии данных одновременно на хосте и в памяти устройства, и по мере необходимости автоматически передаёт обновления между ними.

Кооперативные группы

Ах, как бы хотелось, чтобы можно было написать параллельные алгоритмы для графического процессора единожды — собрать стопку абстракций, обернуть их в шаблоны, а дальше среда выполнения пусть сама разбирается. Но на практике такое удаётся редко. К сожалению, и API CUDA для обращения с кооперативными группами — не исключение.

Он проектировался как единообразная абстракция для координации потоков, не ограниченная пределами одного блока. При этом встроенные функции C++ применяются для назначения сложных алгоритмов GPU, и задействуемая с ними семантика синхронизации достаточно гибкая. Теоретически, так должна решаться серьёзная проблема: пусть все потоки на устройстве время от времени синхронизируются перед продолжением работы — это чрезвычайно важно для итеративных алгоритмов, на которых основаны, например, симуляции физических процессов или решалки.

Просто напомню:

• __syncwarp() обёртывает 32 потока.

• __syncthreads() обёртывает логический блок, включающий 1-1024 потоков.

• На все прочие случаи есть Mastercardкооперативные группы.

Сильнее всего напрашивается такой пример: синхронизируем весь грид в виде многоэтапных итеративных алгоритмов, скажем, при симуляции физического процесса. Для этого Nvidia рекомендует задействовать новую функцию cooperative_groups::sync():

#include <cuda_runtime.h>
#include <cooperative_groups.h>
#include <stdio.h>
#include <math.h>

namespace cg = cooperative_groups;

__device__ float3 compute_force(float3 position_first, float3 position_second) {
    float3 r;
    r.x = position_second.x - position_first.x;
    r.y = position_second.y - position_first.y;
    r.z = position_second.z - position_first.z;

    float squared_distance = r.x * r.x + r.y * r.y + r.z * r.z + 1e-6f; // avoid div by zero
    float reciprocal_distance = rsqrtf(squared_distance);
    float reciprocal_cube = reciprocal_distance * reciprocal_distance * reciprocal_distance;

    constexpr float gravitational_constant = 1.0f;
    float scale = gravitational_constant * reciprocal_cube;
    r.x *= scale;
    r.y *= scale;
    r.z *= scale;
    return r;
}

__global__ void cooperative_kernel(
    float3 *positions_old, float3 *positions_new,
    float3 *velocities_old, float3 *velocities_new,
    size_t count, size_t iterations, float dt) {

    cg::grid_group grid = cg::this_grid();
    size_t idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= count) return;

    for (size_t iter = 0; iter < iterations; ++iter) {
        float3 force = {0.0f, 0.0f, 0.0f};

        // Аккумулируем силы от всех прочих частиц
        for (size_t j = 0; j < count; ++j) {
            if (j == idx) continue;
            float3 f = compute_force(positions_old[idx], positions_old[j]);
            force.x += f.x;
            force.y += f.y;
            force.z += f.z;
        }

        // Обновляем значения скорости и положения
        velocities_new[idx].x = velocities_old[idx].x + force.x * dt;
        velocities_new[idx].y = velocities_old[idx].y + force.y * dt;
        velocities_new[idx].z = velocities_old[idx].z + force.z * dt;
        positions_new[idx].x = positions_old[idx].x + velocities_new[idx].x * dt;
        positions_new[idx].y = positions_old[idx].y + velocities_new[idx].y * dt;
        positions_new[idx].z = positions_old[idx].z + velocities_new[idx].z * dt;

        // Меняем буферы, готовясь к следующей итерации
        grid.sync();
        float3 *temp_pos = positions_old, *temp_vel = velocities_old;
        positions_old = positions_new, positions_new = temp_pos;
        velocities_old = velocities_new, velocities_new = temp_vel;
        grid.sync();
    }
}

int main() {
    cudaError_t err;
    size_t num_particles = 256;
    size_t iterations = 10;
    float dt = 0.01f;
    dim3 block;
    dim3 grid;
    void *kernel_args[7];
    float3 *positions_old = nullptr, *positions_new = nullptr;
    float3 *velocities_old = nullptr, *velocities_new = nullptr;

    // Выделяем память
    err = cudaMallocManaged(&positions_old, num_particles * sizeof(float3));
    if (err != cudaSuccess) goto cleanup;
    err = cudaMallocManaged(&positions_new, num_particles * sizeof(float3));
    if (err != cudaSuccess) goto cleanup;
    err = cudaMallocManaged(&velocities_old, num_particles * sizeof(float3));
    if (err != cudaSuccess) goto cleanup;
    err = cudaMallocManaged(&velocities_new, num_particles * sizeof(float3));
    if (err != cudaSuccess) goto cleanup;

    // Инициализируем значения положения и скорости 
    for (size_t i = 0; i < num_particles; ++i) {
        float theta = (float)i * 0.01f;
        float phi = (float)i * 0.005f;
        float radius = 10.0f + (i % 32) * 0.1f;
        positions_old[i] = {radius * cosf(theta) * sinf(phi), radius * sinf(theta) * sinf(phi), radius * cosf(phi)};
        velocities_old[i] = {0.01f * sinf(phi), 0.01f * cosf(theta), 0.01f * sinf(theta + phi)};
    }

    // Убедимся, что на данном устройстве поддерживается кооперативный запуск потоков
    cudaDeviceProp props;
    cudaGetDeviceProperties(&props, 0);
    if (!props.cooperativeLaunch) {
        fprintf(stderr, "Cooperative launch not supported on this device.\n");
        err = cudaErrorNotSupported;
        goto cleanup;
    }
    cudaStream_t stream;
    err = cudaStreamCreate(&stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to create stream: %s\n", cudaGetErrorString(err));
        goto cleanup;
    }

    block = dim3(256);
    grid = dim3((num_particles + block.x - 1) / block.x);
    kernel_args[0] = &positions_old;
    kernel_args[1] = &positions_new;
    kernel_args[2] = &velocities_old;
    kernel_args[3] = &velocities_new;
    kernel_args[4] = &num_particles;
    kernel_args[5] = &iterations;
    kernel_args[6] = &dt;

    // Запускаем ядро
    err = cudaLaunchCooperativeKernel((void *)cooperative_kernel, grid, block, kernel_args, 0, stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to launch cooperative kernel: %s\n", cudaGetErrorString(err));
        goto cleanup;
    }
    err = cudaStreamSynchronize(stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Kernel execution failed: %s\n", cudaGetErrorString(err));
        goto cleanup;
    }

    // Выводим на экран окончательные положения
    for (size_t i = 0; i < num_particles; ++i)
        printf("Final position[%zu] = (%f, %f, %f)\n", i, positions_old[i].x, positions_old[i].y, positions_old[i].z);

cleanup:
    if (positions_old) cudaFree(positions_old);
    if (positions_new) cudaFree(positions_new);
    if (velocities_old) cudaFree(velocities_old);
    if (velocities_new) cudaFree(velocities_new);
    return (err == cudaSuccess) ? 0 : -1;
}

Обратите внимание, как я заменил cudaLaunchKernel на cudaLaunchCooperativeKernel и добавил к ядру объект cg::grid_group. Если бы мы воспользовались старым некооперативным пусковым API cudaLaunchKernel, то у нас бы получилось:

$ nvcc -o hello_world hello_world.cu && ./hello_world
> Kernel execution failed: unspecified launch failure

Таким образом, следует работать с новым «кооперативным» API, который сразу выглядел многообещающе. При ускоренных связях GPU-GPU и внутриузловых переключателях NVLink можно надеяться, что у нас будут более надёжные примитивы синхронизации для систем со множеством  GPU. В какой-то момент показалось, что это будущее почти наступило: в среде выполнения CUDA ввели cudaLaunchCooperativeKernelMultiDevice и абстракцию cg::multi_grid_group — те самые недостающие звенья, без которых было сложно координировать работу ядер в масштабе множества GPU. Но в CUDA 11.3 оба эти нововведения были признаны устаревшими и таким образом стали одними из самых мимолётных API в истории CUDA.

При всей концептуальной привлекательности кооперативных групп, не думаю, что они смогут внятно масштабироваться. Работа, которой я занимаюсь, в основном происходит на уровне обёртки с использованием __syncwarp(), а если мне требуется подняться на уровень выше, я предпочитаю прибегать к встраиваемому PTX-ассемблеру. Не составляет труда проверить, во что компилируется функция cooperative_groups::sync(): просто используем NVCC с флагом -ptx:

$ nvcc -arch=sm_80 -ptx -o hello_world.ptx hello_world.cu
$ grep -A 1 "barrier.sync" hello_world.ptx

Наши подозрения оправдываются: под капотом это обычная инструкция barrier.sync. Таким образом, если вы хорошо владеете встраиваемым PTX, то можете воспроизвести вышеописанное поведение через заголовок <cooperative_groups.h>. А если вы пишете код, в котором требуется обеспечить максимальную производительность, то он получится не только чище, но и прозрачнее, а также будет лучше поддаваться отладке.

#include <cuda_runtime.h>
#include <stdio.h>
#include <math.h>

__device__ inline void grid_sync_ptx() { asm volatile("barrier.sync 0;" ::); }

__device__ float3 compute_force(float3 position_first, float3 position_second) {
    float3 r;
    r.x = position_second.x - position_first.x;
    r.y = position_second.y - position_first.y;
    r.z = position_second.z - position_first.z;

    float squared_distance = r.x * r.x + r.y * r.y + r.z * r.z + 1e-6f; // избегаем деления на ноль 
    float reciprocal_distance = rsqrtf(squared_distance);
    float reciprocal_cube = reciprocal_distance * reciprocal_distance * reciprocal_distance;

    constexpr float gravitational_constant = 1.0f;
    float scale = gravitational_constant * reciprocal_cube;
    r.x *= scale;
    r.y *= scale;
    r.z *= scale;
    return r;
}

__global__ void cooperative_kernel(float3 *positions_old, float3 *positions_new, float3 *velocities_old, float3 *velocities_new, size_t count,
                                   size_t iterations, float dt) {
    size_t idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= count) return;

    for (size_t iter = 0; iter < iterations; ++iter) {
        float3 force = {0.0f, 0.0f, 0.0f};

        // Аккумулируем силы от всех прочих частиц
        for (size_t j = 0; j < count; ++j) {
            if (j == idx) continue;
            float3 f = compute_force(positions_old[idx], positions_old[j]);
            force.x += f.x;
            force.y += f.y;
            force.z += f.z;
        }

        // Обновляем значения скорости и положения
        velocities_new[idx].x = velocities_old[idx].x + force.x * dt;
        velocities_new[idx].y = velocities_old[idx].y + force.y * dt;
        velocities_new[idx].z = velocities_old[idx].z + force.z * dt;

        positions_new[idx].x = positions_old[idx].x + velocities_new[idx].x * dt;
        positions_new[idx].y = positions_old[idx].y + velocities_new[idx].y * dt;
        positions_new[idx].z = positions_old[idx].z + velocities_new[idx].z * dt;

        grid_sync_ptx();

        // Меняем буферы, готовясь к следующей итерации
        float3 *temp_pos = positions_old, *temp_vel = velocities_old;
        positions_old = positions_new, positions_new = temp_pos;
        velocities_old = velocities_new, velocities_new = temp_vel;

        grid_sync_ptx();
    }
}

int main() {
    cudaError_t err;
    size_t num_particles = 256;
    size_t iterations = 10;
    float dt = 0.01f;

    float3 *positions_old = nullptr, *positions_new = nullptr;
    float3 *velocities_old = nullptr, *velocities_new = nullptr;

    dim3 block;
    dim3 grid;
    void *kernel_args[7];

    // Выделяем память
    err = cudaMallocManaged(&positions_old, num_particles * sizeof(float3));
    if (err != cudaSuccess) goto cleanup;
    err = cudaMallocManaged(&positions_new, num_particles * sizeof(float3));
    if (err != cudaSuccess) goto cleanup;
    err = cudaMallocManaged(&velocities_old, num_particles * sizeof(float3));
    if (err != cudaSuccess) goto cleanup;
    err = cudaMallocManaged(&velocities_new, num_particles * sizeof(float3));
    if (err != cudaSuccess) goto cleanup;

    for (size_t i = 0; i < num_particles; ++i) {
        float theta = (float)i * 0.01f;
        float phi = (float)i * 0.005f;
        float radius = 10.0f + (i % 32) * 0.1f;
        positions_old[i] = {radius * cosf(theta) * sinf(phi), radius * sinf(theta) * sinf(phi), radius * cosf(phi)};
        velocities_old[i] = {0.01f * sinf(phi), 0.01f * cosf(theta), 0.01f * sinf(theta + phi)};
    }

    cudaDeviceProp props;
    cudaGetDeviceProperties(&props, 0);
    if (!props.cooperativeLaunch) {
        fprintf(stderr, "Cooperative launch not supported on this device.\n");
        err = cudaErrorNotSupported;
        goto cleanup;
    }

    cudaStream_t stream;
    err = cudaStreamCreate(&stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to create stream: %s\n", cudaGetErrorString(err));
        goto cleanup;
    }

    block = dim3(256);
    grid = dim3((num_particles + block.x - 1) / block.x);
    kernel_args[0] = &positions_old;
    kernel_args[1] = &positions_new;
    kernel_args[2] = &velocities_old;
    kernel_args[3] = &velocities_new;
    kernel_args[4] = &num_particles;
    kernel_args[5] = &iterations;
    kernel_args[6] = &dt;

    err = cudaLaunchKernel((void *)cooperative_kernel, grid, block, kernel_args, 0, stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Failed to launch cooperative kernel: %s\n", cudaGetErrorString(err));
        goto cleanup;
    }

    err = cudaStreamSynchronize(stream);
    if (err != cudaSuccess) {
        fprintf(stderr, "Kernel execution failed: %s\n", cudaGetErrorString(err));
        goto cleanup;
    }

    for (size_t i = 0; i < num_particles; ++i)
        printf("Final position[%zu] = (%f, %f, %f)\n", i, positions_old[i].x, positions_old[i].y, positions_old[i].z);

cleanup:
    if (positions_old) cudaFree(positions_old);
    if (positions_new) cudaFree(positions_new);
    if (velocities_old) cudaFree(velocities_old);
    if (velocities_new) cudaFree(velocities_new);
    return (err == cudaSuccess) ? 0 : -1;
}

Кстати, пока никто не видел, мы запустились по старинке при помощи  cudaLaunchKernel, а не с cudaLaunchCooperativeKernel. Среда выполнения на это жаловаться не будет, при условии, что у вас на устройстве поддерживается barrier.sync. Если вы азартны, то PTX найдётся для вас и целый спектр других барьеров.

API драйверов CUDA

Наконец, незаслуженно обходят вниманием ещё более низкоуровневый API драйверов CUDA, который может быть исключительно полезен на практике. Да, кода в нём многовато, зато он предоставляет вам полный контроль над загрузкой и запуском ядра, в том числе, поддерживает динамическую загрузку PTX, CUBIN или SASS во время выполнения. Немного забегая вперёд, рекомендую полностью отграничить код ядра от кода хоста, причём, использовать для них два отдельных компилятора и обустроить между ними стабильный ABI. Вот как упрощённом виде мог бы выглядеть наш код хоста на C99:

#include <cuda.h>
#include <stdio.h>
#include <math.h>

#define CUDA_CHECK(err)                             \
    if (err != CUDA_SUCCESS) {                      \
        const char *msg;                            \
        cuGetErrorString(err, &msg);                \
        fprintf(stderr, "CUDA error: %s\n", msg);   \
        goto cleanup;                               \
    }

int main() {
    CUresult err;
    size_t num_particles = 256;
    size_t iterations = 10;
    float dt = 0.01f;

    CUdevice device;
    CUcontext context = NULL;
    CUmodule module = NULL;
    CUfunction kernel;
    CUstream stream = NULL;
    float *positions_old = NULL, *positions_new = NULL;
    float *velocities_old = NULL, *velocities_new = NULL;
    void *kernel_args[7];

    // Инициализация CUDA
    err = cuInit(0);
    CUDA_CHECK(err);
    err = cuDeviceGet(&device, 0);
    CUDA_CHECK(err);
    err = cuCtxCreate(&context, 0, device);
    CUDA_CHECK(err);
    err = cuStreamCreate(&stream, CU_STREAM_DEFAULT);
    CUDA_CHECK(err);

    // Загрузка PTX "байт-кода", который впоследствии будет динамически компилироваться в SASS
    err = cuModuleLoad(&module, "hello_world.ptx");
    CUDA_CHECK(err);
    err = cuModuleGetFunction(&kernel, module, "cooperative_kernel");
    CUDA_CHECK(err);

    // Выделение управляемой памяти для положений и скоростей 
    size_t buffer_size = num_particles * sizeof(float) * 3;
    err = cuMemAllocManaged((CUdeviceptr *)&positions_old, buffer_size, CU_MEM_ATTACH_GLOBAL);
    CUDA_CHECK(err);
    err = cuMemAllocManaged((CUdeviceptr *)&positions_new, buffer_size, CU_MEM_ATTACH_GLOBAL);
    CUDA_CHECK(err);
    err = cuMemAllocManaged((CUdeviceptr *)&velocities_old, buffer_size, CU_MEM_ATTACH_GLOBAL);
    CUDA_CHECK(err);
    err = cuMemAllocManaged((CUdeviceptr *)&velocities_new, buffer_size, CU_MEM_ATTACH_GLOBAL);
    CUDA_CHECK(err);

    // Инициализация положений и скоростей 
    for (size_t i = 0; i < num_particles; ++i) {
        float theta = (float)i * 0.01f;
        float phi = (float)i * 0.005f;
        float radius = 10.0f + (i % 32) * 0.1f;
        positions_old[3 * i + 0] = radius * cosf(theta) * sinf(phi);
        positions_old[3 * i + 1] = radius * sinf(theta) * sinf(phi);
        positions_old[3 * i + 2] = radius * cosf(phi);
        velocities_old[3 * i + 0] = 0.01f * sinf(phi);
        velocities_old[3 * i + 1] = 0.01f * cosf(theta);
        velocities_old[3 * i + 2] = 0.01f * sinf(theta + phi);
    }

    kernel_args[0] = &positions_old;
    kernel_args[1] = &positions_new;
    kernel_args[2] = &velocities_old;
    kernel_args[3] = &velocities_new;
    kernel_args[4] = &num_particles;
    kernel_args[5] = &iterations;
    kernel_args[6] = &dt;

    // Запуск ядра
    int threads_per_block = 256;
    int blocks_per_grid = (num_particles + threads_per_block - 1) / threads_per_block;
    err = cuLaunchKernel(kernel,
                         blocks_per_grid, 1, 1,
                         threads_per_block, 1, 1,
                         0, stream,
                         kernel_args, NULL);
    CUDA_CHECK(err);
    err = cuStreamSynchronize(stream);
    CUDA_CHECK(err);

    // Логирование окончательных положений
    for (size_t i = 0; i < num_particles; ++i)
        printf("Final position[%zu] = (%f, %f, %f)\n", i,
               positions_old[3 * i + 0],
               positions_old[3 * i + 1],
               positions_old[3 * i + 2]);

cleanup:
    if (stream) cuStreamDestroy(stream);
    if (positions_old) cuMemFree((CUdeviceptr)positions_old);
    if (positions_new) cuMemFree((CUdeviceptr)positions_new);
    if (velocities_old) cuMemFree((CUdeviceptr)velocities_old);
    if (velocities_new) cuMemFree((CUdeviceptr)velocities_new);
    if (module) cuModuleUnload(module);
    if (context) cuCtxDestroy(context);
    return (err == CUDA_SUCCESS) ? 0 : -1;
}

В данном случае особого внимания требует лишь декорирование имён. Чтобы этот механизм работал, убедитесь, что объявление ядра у вас в файле .cu обёрнуто в extern "C":

extern "C" __global__ void cooperative_kernel(...);

Уладив это, можно компилировать код GPU в PTX при помощи NVCC, а код хоста при помощи GCC — совершенно независимо:

Заключение

Конечно, у нас получилось гораздо больше кода, чем в исходном примере с <<<1, 1>>> — но обычно так и бывает, если стараешься всё делать правильно. При этом отмечу, что в последние годы появилось множество инструментов, упрощающих прототипирование кода CUDA, в том числе, разнообразные DSL и компиляторы, которыми NVIDIA любит козырнуть. Но, если готовить и выпускать ядра для продакшена так, как описано здесь, то они получаются удивительно стабильными и работают по 10 и более лет.

Изменилась лишь сложность самих ядер — теперь в них меньше распараллеливаются данные, и ядра сближаются по принципу работы с конкурентными алгоритмами для ЦП. Здесь применяются и атомарные операции, и редукция на уровне варпов, и во всё это запечена логика тензорных ядер.  Эти реализации выглядят и ощущаются по-разному в каждом новом поколении GPU, в чём мне довелось убедиться на собственной шкуре при портировании других библиотек.