Первые ядра: от сложения чисел к векторам

1. Разминка: сложить два числа

Задача абсурдно мала для видеокарты — тем и хороша: в программе не останется ничего, кроме протокола из C03, и каждый шаг видно под лупой.

#include <cstdio>

__global__ void add(int* a, int* b, int* c)
{
    *c = *a + *b;                 // один тред складывает два числа
}

int main()
{
    int a = 5, b = 10, c = 0;     // хост-копии
    int *da, *db, *dc;            // device-указатели
    int size = sizeof(int);

    cudaMalloc((void**)&da, size);            // 1. выделить на устройстве
    cudaMalloc((void**)&db, size);
    cudaMalloc((void**)&dc, size);

    cudaMemcpy(da, &a, size, cudaMemcpyHostToDevice);   // 2. вход туда
    cudaMemcpy(db, &b, size, cudaMemcpyHostToDevice);

    add<<<1, 1>>>(da, db, dc);                          // 3. посчитать

    cudaMemcpy(&c, dc, size, cudaMemcpyDeviceToHost);   // 4. результат обратно

    printf("%d + %d = %d\n", a, b, c);                 // 5 + 10 = 15
    cudaFree(da); cudaFree(db); cudaFree(dc);           // 5. освободить
    return 0;
}

Обратите внимание: между запуском ядра и cudaMemcpy нет cudaDeviceSynchronize, хотя запуск асинхронный. Это не ошибка: cudaMemcpy сам ждёт завершения всей предыдущей работы устройства — синхронизация встроена в него. Явная синхронизация нужна, когда после ядра идёт не копирование, а, например, печать с хоста (как в C01) или замер времени (C06).

Ядро выполнил один тред одного блока. Видеокарта с тысячами ядер простаивала — время её загрузить.

2. Сложение векторов: три конфигурации

Теперь складываем поэлементно два массива по n чисел. Работы n штук, и раздать её можно по-разному. У ДВФУ эти варианты идут как три отдельные программы; сведём их в таблицу — меняется только индекс и конфигурация запуска:

Формула третьей строки — глобальный индекс треда — самая используемая строчка в CUDA. Читается она так: пропустить все предыдущие блоки (blockIdx.x штук по blockDim.x тредов) и добавить свой номер в текущем блоке. Вы уже щёлкали по ней в тренажёре C02; теперь она встала на рабочее место:

__global__ void vecAdd(const float* a, const float* b, float* c, int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n)                    // защита: лишние треды не лезут за край массива
        c[i] = a[i] + b[i];
}

Откуда лишние треды? Из арифметики. Пусть n = 1000, а блок — 256 тредов. Блоков нужно 1000 / 256 = 3.9, а блоки бывают только целыми. Три блока дадут 768 тредов — мало; четыре — 1024, на 24 больше, чем элементов. Поэтому блоков берут с округлением вверх, а хвостовые треды отсекает то самое if (i < n):

const int block = 256;
int grid = (n + block - 1) / block;   // целочисленное «округление вверх»:
                                       // (1000 + 255) / 256 = 1255/256 = 4
vecAdd<<<grid, block>>>(da, db, dc, n);

Трюк (n + block − 1) / block стоит понять один раз навсегда: прибавка block − 1 «доталкивает» любой неполный остаток до следующего целого, а полные значения не трогает (при n = 1024: (1024 + 255)/256 = 4 — ровно). Полная программа вокруг ядра — тот же протокол из раздела 1, только size = n · sizeof(float):

#include <cstdio>

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

int main()
{
    const int n = 1000;
    int size = n * sizeof(float);

    // хост: обычные массивы, заполняем тестовыми значениями
    float* a = new float[n];
    float* b = new float[n];
    float* c = new float[n];
    for (int i = 0; i < n; i++) { a[i] = i; b[i] = 2 * i; }

    // устройство: выделить и завезти вход
    float *da, *db, *dc;
    cudaMalloc((void**)&da, size);
    cudaMalloc((void**)&db, size);
    cudaMalloc((void**)&dc, size);
    cudaMemcpy(da, a, size, cudaMemcpyHostToDevice);
    cudaMemcpy(db, b, size, cudaMemcpyHostToDevice);

    // запуск: блоков с округлением вверх
    const int block = 256;
    int grid = (n + block - 1) / block;
    vecAdd<<<grid, block>>>(da, db, dc, n);

    // результат и проверка: c[i] должно равняться 3i
    cudaMemcpy(c, dc, size, cudaMemcpyDeviceToHost);
    bool ok = true;
    for (int i = 0; i < n; i++)
        if (c[i] != 3.0f * i) { printf("ошибка в %d: %f\n", i, c[i]); ok = false; break; }
    if (ok) printf("все %d элементов верны\n", n);

    cudaFree(da); cudaFree(db); cudaFree(dc);
    delete[] a; delete[] b; delete[] c;
    return 0;
}

Проверка результата на хосте — не паранойя, а привычка: ядро, вычислившее мусор, не сообщает об этом никак. Дешёвый цикл сравнения после каждого нового ядра экономит часы поиска «куда делись данные».

3. Когда цикл можно раздать тредам

Сложение векторов легло на GPU без сопротивления. Почему? Взгляните на его последовательную версию:

for (int i = 0; i < n; i++)
    c[i] = a[i] + b[i];          // итерация i не трогает ничьих чужих данных

Каждая итерация читает только a[i] и b[i], пишет только в c[i] — и ей безразлично, выполнены ли остальные итерации, идут ли они сейчас или ещё не начинались. Это критерий независимости итераций, и он же — критерий «цикл можно раздать тредам»: строка тела цикла превращается в ядро, счётчик i — в глобальный индекс, и порядок выполнения перестаёт иметь значение (а C02 научил: порядка и не будет). Такой паттерн называют map: одна независимая операция на элемент.

Теперь антипример — накопленная сумма:

for (int i = 1; i < n; i++)
    x[i] = x[i - 1] + a[i];     // итерация i ЖДЁТ результат итерации i-1

Итерация i читает x[i−1] — то, что записала предыдущая итерация. Раздать этот цикл тредам «в лоб» нельзя: тред i запустится когда угодно и почти наверняка прочитает x[i−1] до того, как тред i−1 туда что-то запишет. Получится мусор — причём каждый запуск разный. Это рекуррентность, перенос зависимости между итерациями.

Между чистым map и безнадёжной рекуррентностью лежат два важных промежуточных паттерна. Reduce — свёртка массива в одно число (сумма, максимум): итерации зависят через общий аккумулятор, но операция ассоциативна, и зависимость можно перестроить в дерево частичных сумм — этим займётся C05. Scan — префиксные суммы: выглядит как рекуррентность, но тоже параллелится специальным алгоритмом. Полная классификация с примерами на Python и разбором того, что умеет и чего не умеет автоматический параллелизатор, — в конспекте P02 соседней линии; для линии C достаточно правила:

4. Тренажёр: индексация до автоматизма

Контрольные вопросы

• Выведите формулу глобального индекса и объясните каждый множитель.
  i = blockIdx.x · blockDim.x + threadIdx.x. Первое слагаемое пропускает все предыдущие блоки: их blockIdx.x штук, в каждом blockDim.x тредов. Второе — номер треда внутри своего блока. Вместе получается сквозная нумерация тредов грида, которой удобно нумеровать элементы данных.
• Зачем нужна защита if (i < n) и когда её можно опустить?
  Число тредов кратно размеру блока и обычно больше n: хвостовые треды получают индексы за границей массива. Защита не даёт им читать и писать чужую память. Опустить можно, только если n гарантированно кратно размеру блока — но привычнее и безопаснее писать всегда.
• Почему (n + block − 1) / block даёт округление вверх?
  Целочисленное деление отбрасывает остаток. Прибавка block − 1 переводит любое n с ненулевым остатком через следующую границу кратности (остаток r ≥ 1 превращается в r + block − 1 ≥ block), а кратные n не сдвигает: их остаток 0, и добавка не дотягивает до следующего целого.
• Чем плохи конфигурации «n блоков по 1 треду» и «1 блок по n тредов»?
  Блок из одного треда занимает целый варп: 31 из 32 мест простаивает. Один блок на всё — упирается в лимит 1024 треда и работает на единственном SM, остальные мультипроцессоры бездельничают. Рабочая схема — много блоков по 128–512 тредов: масштабируется и на длину данных, и на число SM карты.
• Сформулируйте критерий «цикл можно раздать тредам» и приведите контрпример.
  Итерации независимы: каждая читает и пишет только свои данные и не использует результатов других итераций — тогда порядок выполнения безразличен, и тела итераций можно выполнять одновременно (паттерн map). Контрпример — рекуррентность x[i] = x[i−1] + a[i]: итерация ждёт результата предыдущей, наивная параллелизация даёт мусор, разный от запуска к запуску.
• Почему между запуском ядра и cudaMemcpy не нужен cudaDeviceSynchronize?
  cudaMemcpy сам синхронизируется: он не начнёт копирование, пока устройство не закончит предыдущую работу, и не вернёт управление, пока не скопирует. Явная синхронизация нужна, когда после ядра идёт не CUDA-вызов, а хостовый код — печать, замер времени системными часами.

Источники

1. An Even Easier Introduction to CUDA // NVIDIA Developer Blog : [сайт]. — URL: https://developer.nvidia.com/blog/even-easier-introduction-cuda/ (дата обращения: 09.07.2026).
2. CUDA C++ Programming Guide. Kernels; Thread Hierarchy // NVIDIA Docs : [сайт]. — URL: https://docs.nvidia.com/cuda/cuda-c-programming-guide/ (дата обращения: 09.07.2026).