Доступ к памяти делает ядро ​​CUDA чрезвычайно медленным

Сначала несколько замечаний по вашим фактическим вычислениям:

• Вы сравните sqrtf(x) < sqrtf(3). Корни дорогие. Просто сравните x < 3.f
• Даже если вы хотите сохранить квадратный корень, чтобы избежать переполнения диапазона с плавающей запятой (вероятно, это не проблема), не используйте sqrt(pow(x, 2)+...), в этом отношении не используйте pow только для возведения в квадрат. Используйте hypotf для 2D или norm3df для 3D векторов.
• Вы сохраняете последнее значение в вашей палитре, которое ниже предела. Очень похоже, что вы хотите подобрать лучший цвет

Теперь давайте проанализируем ваши обращения к памяти:

индекс фрагмента

Давайте посмотрим на fragment_idx = 3*gridDim.y*blockIdx.x+3*blockIdx.y: вы не принимаете во внимание threadIdx.x и threadIdx.y. Это ваша основная проблема: многие потоки действуют на одни и те же входные и выходные данные. Вы, вероятно, хотите этого: fragment_idx = 3 * (threadIdx.y * blockDim.x + threadIdx.x)

вход

Итак, вы загружаете 3 поплавка. Во-первых, почему вы перезагружаете его внутри цикла, если он не зависит от итерации цикла? Я предполагаю, что компилятор спасает вас от этого доступа, но не привыкайте делать это.

Во-вторых, ваш шаблон доступа не объединен должным образом, поскольку а) это 3 независимых доступа и б) CUDA не может объединить доступы к float3 векторам, даже если вы сделали это правильно. Прочтите раздел 9.2.1 Объединенный доступ к глобальной памяти Руководства по передовому опыту. Для повышения производительности у вас есть два варианта:

1. Вы добавляете 1 число с плавающей запятой на fragment_idx, чтобы вы могли загрузить все это как float4
2. Вы переносите свой входной массив из матрицы Nx3 в матрицу 3xN

палитра

Та же проблема с доступом к 3 поплавкам. Кроме того, теперь каждый поток считывает одни и те же значения, поскольку c не зависит от индекса потока. По крайней мере, доступ должен проходить через функцию __ldg для использования кеша L1. Предпочтительно предварительно загрузить палитру в общую память.

выход

Доступ для записи имеет ту же проблему, что и несвязанный доступ. Кроме того, поскольку best_c варьируется в зависимости от потока, доступ для чтения к palette является случайным. Вы должны были загрузить значения palette раньше в свой цикл. Просто сохраните лучшее значение палитры в локальной переменной и повторно используйте его для сохранения вывода в конце.

Методология

Два замечания:

1. Постарайтесь сделать свой код действительным, прежде чем делать его быстрым. Это бы поймало fragment_idx
2. Если вы упростите код, например удалив вывод, компилятор с радостью оптимизирует большую часть вашего кода. Это не то, как вы делаете правильную оценку производительности. Используйте профайлер. CUDA поставляется с очень хорошими

Минимальное исправление

Это самый простой код для устранения проблем. Это не решает проблемы с загрузкой переменных vector3 и не использует разделяемую память. Это требует более серьезных изменений

__device__ float sqr_norm(float3 a, float3 b) {
    a.x -= b.x, a.y -= b.y, a.z -= b.z;
    return a.x * a.x + a.y * a.y + a.z * a.z;
}
__global__ void fragment_shader(int palette_lim,
          const float *palette, const float *input,
          float *output) {
    int fragment_idx = 3 * (threadIdx.y * blockDim.x + threadIdx.x);
    /* TODO: Switch to float4 for better memory access patterns */
    float3 inputcolor = make_float3(
          input[fragment_idx], input[fragment_idx + 1], input[fragment_idx + 2]);
    float min_dist_sqr = 3.f;
    /* The old code always used color index 0 if there was no fit */
    float3 best_color = make_float3(
          __ldg(palette), __ldg(palette + 1), __ldg(palette + 2));
    float best_dist = sqr_norm(best_color, inputcolor);
    for(int c = 1; c < palette_lim; c++) {
        /* TODO: Prefetch into shared memory */
        float3 color = make_float3(
              __ldg(palette + c), __ldg(palette + c + 1), __ldg(palette + c + 2));
        float dist = sqr_norm(color, inputcolor);
        /* Since we always used color 0 in the old code,
         * the min_dist is somewhat pointless */
        if (dist < min_dist_sqr && dist < best_dist) {
            best_color = color;
            best_dist = dist;
        }
    }
    output[fragment_idx] = best_color.x;
    output[fragment_idx + 1] = best_color.y;
    output[fragment_idx + 2] = best_color.z;
}

Обширное исправление

Вот более обширная переработка:

1. Все массивы изменены на float4 (RGBA вместо RGB). Дополнительный канал игнорируется при вычислении расстояния, но он распространяется. Обычно кто-то пытается использовать значение для чего-то, например. вы можете сохранить значение расстояния там
2. Общая память используется для буферизации цветовой палитры. Требования к динамической общей памяти изложены в комментариях к коду.

__device__ float sqr_dist_rgb(float4 a, float4 b) {
    a.x -= b.x, a.y -= b.y, a.z -= b.z;
    return a.x * a.x + a.y * a.y + a.z * a.z;
}
__global__ void fragment_shader(int palette_lim,
          const float4 *palette, const float4 *input,
          float4 *output) {
    /* Call with dynamic shared memory:
     * 2 * sizeof(float4) * blockDim.x * blockDim.y */
    extern __shared__ float4 colorbuf[];
    const int buf_size = blockDim.x * blockDim.y;
    const int buf_idx = threadIdx.y * blockDim.x + threadIdx.x;
    const int fragment_idx = threadIdx.y * blockDim.x + threadIdx.x;
    const float4 inputcolor = input[fragment_idx];
    float4 best_color =  __ldg(palette);
    const float min_dist_sqr = 3.f;
    float best_dist = sqr_dist_rgb(best_color, inputcolor);
    for(int cb = 0, b = 0; cb < palette_lim; b ^= 1, cb += buf_size) {
        /* We use a double buffer scheme to reduce the __syncthreads calls */
        float4* cur_buf = b ? colorbuf + buf_size : colorbuf;
        if (cb + buf_idx < palette_lim)
            cur_buf[buf_idx] = __ldg(palette + cb + buf_idx);
        __syncthreads();
        const int n = min(buf_size, palette_lim - cb);
        for(int c = 0; c < n; c++) {
            float4 color = cur_buf[c];
            float dist = sqr_dist_rgb(color, inputcolor);
            if (dist < min_dist_sqr && dist < best_dist) {
                best_color = color;
                best_dist = dist;
            }
        }
    }
    output[fragment_idx] = best_color;
}

Алгоритмические улучшения

Для больших палитр такой поиск методом грубой силы неоптимален. Алгоритмы пространственного индекса могут делать то же самое, но быстрее. Классической структурой для этого будет KD-дерево. Если вы будете искать это, вы найдете некоторые реализации CUDA, которые, возможно, стоит проверить.