Как оптимизировать сито CUDA Eratosthenes

Question

Я новичок в CUDA. Чтобы мои руки были грязными, я попробовал написать сито Эратосфена (для нахождения всех простых чисел до некоторого числа n).

Есть несколько вещей, которые я должен был сделать, чтобы заставить его работать, что, похоже, не обязательно было. Мне любопытно, знает ли кто-нибудь более естественный (и все же CUDA-оптимизированный) подход.

1. Чтобы взять записи, помеченные как простые в массиве isPrime, мне пришлось выполнять два отдельных вызова ядра. Первый подсчитывает количество простых чисел в каждом потоковом блоке и присваивает каждой записи i число простых чисел в этом блоке меньше i. Затем я должен сделать второй звонок, чтобы добавить число простых чисел во всех предыдущих блоках, чтобы получить окончательный индекс.
2. Но это даже хуже, потому что, чтобы избежать кучи параллельных чтений, мне пришлось хранить количество простых чисел в блоке в отдельном массиве по каждому из индексов THREADS_PER_BLOCK, эффективно удваивая требуемую память для алгоритма. Кажется, должен быть способ, чтобы все потоки читали одно и то же значение для каждого блока, а не копировали его столько раз.
3. Несмотря на все это, проблема одновременного чтения в методе clearMultiples. Особенно для небольших простых чисел, таких как 2 и 3, каждая нить должна читать значение. Не существует ли способа справиться с этим?

Может ли кто-нибудь взглянуть на мой код и рассказать мне, есть ли что-то очевидное, что я мог бы сделать, это было бы проще или эффективнее?

Есть ли что-то, что я делаю, это особенно неэффективно (кроме того, распечатывать все штрихи в конце курса)?

Нужно ли вызывать синхронизацию после каждого вызова ядра?

Нужно ли мне синхронизировать и после memcpy?

И наконец, как я устанавливаю THREADS_PER_BLOCK на 512, это не работает?

Спасибо

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

#define MAX_BLOCKS 256 
#define THREADS_PER_BLOCK 256 //Must be a power of 2 
#define BLOCK_SPACE 2 * THREADS_PER_BLOCK 

__global__ void initialize(int* isPrime, int n) { 
    int idx = blockIdx.x * THREADS_PER_BLOCK + threadIdx.x; 
    int step = gridDim.x * THREADS_PER_BLOCK; 
    int i; 
    for (i = idx; i <= 1; i += step) { 
     isPrime[i] = 0; 
    } 
    for (; i < n; i += step) { 
     isPrime[i] = 1; 
    } 
} 

__global__ void clearMultiples(int* isPrime, int* primeList, int startInd, 
     int endInd, int n) { 
    int yidx = blockIdx.y * blockDim.y + threadIdx.y; 
    int xidx = blockIdx.x * blockDim.x + threadIdx.x; 
    int ystep = gridDim.y * blockDim.y; 
    int xstep = gridDim.x * blockDim.x; 
    for (int pnum = startInd + yidx; pnum < endInd; pnum += ystep) { 
     int p = primeList[pnum]; 
     int pstart = p * (p + xidx); 
     int pstep = p * xstep; 
     for (int i = pstart; i < n; i += pstep) { 
      isPrime[i] = 0; 
     } 
    } 
} 

__device__ void makeCounts(int* isPrime, int* addend, int start, int stop) { 
    __shared__ int tmpCounts[BLOCK_SPACE]; 
    __shared__ int dumbCounts[BLOCK_SPACE]; 
    int idx = threadIdx.x; 
    tmpCounts[idx] = ((start + idx) < stop) ? isPrime[start + idx] : 0; 
    __syncthreads(); 
    int numEntries = THREADS_PER_BLOCK; 
    int cstart = 0; 
    while (numEntries > 1) { 
     int prevStart = cstart; 
     cstart += numEntries; 
     numEntries /= 2; 
     if (idx < numEntries) { 
      int i1 = idx * 2 + prevStart; 
      tmpCounts[idx + cstart] = tmpCounts[i1] + tmpCounts[i1 + 1]; 
     } 
     __syncthreads(); 
    } 
    if (idx == 0) { 
     dumbCounts[cstart] = tmpCounts[cstart]; 
     tmpCounts[cstart] = 0; 
    } 
    while (cstart > 0) { 
     int prevStart = cstart; 
     cstart -= numEntries * 2; 
     if (idx < numEntries) { 
      int v1 = tmpCounts[idx + prevStart]; 
      int i1 = idx * 2 + cstart; 
      tmpCounts[i1 + 1] = tmpCounts[i1] + v1; 
      tmpCounts[i1] = v1; 
      dumbCounts[i1] = dumbCounts[i1 + 1] = dumbCounts[idx + prevStart]; 
     } 
     numEntries *= 2; 
     __syncthreads(); 
    } 
    if (start + idx < stop) { 
     isPrime[start + idx] = tmpCounts[idx]; 
     addend[start + idx] = dumbCounts[idx]; 
    } 
} 

__global__ void createCounts(int* isPrime, int* addend, int lb, int ub) { 
    int step = gridDim.x * THREADS_PER_BLOCK; 
    for (int i = lb + blockIdx.x * THREADS_PER_BLOCK; i < ub; i += step) { 
     int start = i; 
     int stop = min(i + step, ub); 
     makeCounts(isPrime, addend, start, stop); 
    } 
} 

__global__ void sumCounts(int* isPrime, int* addend, int lb, int ub, 
     int* totalsum) { 
    int idx = blockIdx.x; 
    int s = 0; 
    for (int i = lb + idx; i < ub; i += THREADS_PER_BLOCK) { 
     isPrime[i] += s; 
     s += addend[i]; 
    } 
    if (idx == 0) { 
     *totalsum = s; 
    } 
} 

__global__ void condensePrimes(int* isPrime, int* primeList, int lb, int ub, 
     int primeStartInd, int primeCount) { 
    int idx = blockIdx.x * THREADS_PER_BLOCK + threadIdx.x; 
    int step = gridDim.x * THREADS_PER_BLOCK; 
    for (int i = lb + idx; i < ub; i += step) { 
     int term = isPrime[i]; 
     int nextTerm = i + 1 == ub ? primeCount : isPrime[i + 1]; 
     if (term < nextTerm) { 
      primeList[primeStartInd + term] = i; 
     } 
    } 
} 

int main(void) { 
    printf("Enter upper bound:\n"); 
    int n; 
    scanf("%d", &n); 

    int *isPrime, *addend, *numPrimes, *primeList; 
    cudaError_t t = cudaMalloc((void**) &isPrime, n * sizeof(int)); 
    assert(t == cudaSuccess); 
    t = cudaMalloc(&addend, n * sizeof(int)); 
    assert(t == cudaSuccess); 
    t = cudaMalloc(&numPrimes, sizeof(int)); 
    assert(t == cudaSuccess); 
    int primeBound = 2 * n/log(n); 
    t = cudaMalloc(&primeList, primeBound * sizeof(int)); 
    assert(t == cudaSuccess); 
    int numBlocks = min(MAX_BLOCKS, 
      (n + THREADS_PER_BLOCK - 1)/THREADS_PER_BLOCK); 
    initialize<<<numBlocks, THREADS_PER_BLOCK>>>(isPrime, n); 
    t = cudaDeviceSynchronize(); 
    assert(t == cudaSuccess); 

    int bound = (int) ceil(sqrt(n)); 

    int lb; 
    int ub = 2; 
    int primeStartInd = 0; 
    int primeEndInd = 0; 

    while (ub < n) { 
     if (primeEndInd > primeStartInd) { 
      int lowprime; 
      t = cudaMemcpy(&lowprime, primeList + primeStartInd, sizeof(int), 
        cudaMemcpyDeviceToHost); 
      assert(t == cudaSuccess); 
      int numcols = n/lowprime; 
      int numrows = primeEndInd - primeStartInd; 
      int threadx = min(numcols, THREADS_PER_BLOCK); 
      int thready = min(numrows, THREADS_PER_BLOCK/threadx); 
      int blockx = min(numcols/threadx, MAX_BLOCKS); 
      int blocky = min(numrows/thready, MAX_BLOCKS/blockx); 
      dim3 gridsize(blockx, blocky); 
      dim3 blocksize(threadx, thready); 
      clearMultiples<<<gridsize, blocksize>>>(isPrime, primeList, 
        primeStartInd, primeEndInd, n); 
      t = cudaDeviceSynchronize(); 
      assert(t == cudaSuccess); 
     } 
     lb = ub; 
     ub *= 2; 
     if (lb >= bound) { 
      ub = n; 
     } 
     numBlocks = min(MAX_BLOCKS, 
       (ub - lb + THREADS_PER_BLOCK - 1)/THREADS_PER_BLOCK); 

     createCounts<<<numBlocks, THREADS_PER_BLOCK>>>(isPrime, addend, lb, ub); 
     t = cudaDeviceSynchronize(); 
     assert(t == cudaSuccess); 

     sumCounts<<<THREADS_PER_BLOCK, 1>>>(isPrime, addend, lb, ub, numPrimes); 
     t = cudaDeviceSynchronize(); 
     assert(t == cudaSuccess); 

     int primeCount; 
     t = cudaMemcpy(&primeCount, numPrimes, sizeof(int), 
       cudaMemcpyDeviceToHost); 
     assert(t == cudaSuccess); 
     assert(primeCount > 0); 

     primeStartInd = primeEndInd; 
     primeEndInd += primeCount; 
     condensePrimes<<<numBlocks, THREADS_PER_BLOCK>>>(isPrime, primeList, lb, 
       ub, primeStartInd, primeCount); 
     t = cudaDeviceSynchronize(); 
     assert(t == cudaSuccess); 
    } 
    int finalprimes[primeEndInd]; 
    t = cudaMemcpy(finalprimes, primeList, primeEndInd * sizeof(int), 
      cudaMemcpyDeviceToHost); 
    assert(t == cudaSuccess); 

    t = cudaFree(isPrime); 
    assert(t == cudaSuccess); 
    t = cudaFree(addend); 
    assert(t == cudaSuccess); 
    t = cudaFree(numPrimes); 
    assert(t == cudaSuccess); 
    t = cudaFree(primeList); 
    assert(t == cudaSuccess); 
    for (int i = 0; i < primeEndInd; i++) { 
     if (i % 16 == 0) 
      printf("\n"); 
     else 
      printf(" "); 
     printf("%4d", finalprimes[i]); 
    } 
    printf("\n"); 
    return 0; 
} 

Answer 1

Отвечая на некоторые ваши вопросы.

1. Исправьте свою проверку ошибок, как определено в комментариях.

2. определите, что вы подразумеваете под «одновременным чтением». Вы обеспокоены этим, но я не уверен, что вы подразумеваете под этим.

Нужно вызвать синхронизацию после каждого вызова ядра?

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

Нужно ли мне синхронизировать и после memcpy?

Нет, cudaMemcpy синхронный характер (это заставит все вызовы Cuda в том же потоке, чтобы завершить до его начала, и он не будет возвращать управление принимающей нить, пока не будет завершена.) Если вы не Не нужно блокировать характеристику (не возвращая управление потоку хоста до завершения), тогда вы можете использовать версию вызова cudaMemcpyAsync. Вы должны использовать потоки, чтобы обойти поведение принудительного завершения всех предыдущих вызовов cuda.

И наконец, как я установил THREADS_PER_BLOCK на 512, это не сработает?

Пожалуйста, определите, что вы подразумеваете под «это не работает». Я скомпилировал ваш код с THREADS_PER_BLOCK 512 и 256, а для верхней границы 1000 он дал одинаковый результат в каждом случае.