Оптимизация N-Body CUDA

Question

Я разрабатываю алгоритм N-body в CUDA и хотел бы выучить некоторые советы и рекомендации по оптимизации.

Мне удалось заставить 16384 тела работать на 20Flops вNVIDIA Geforce GTX 260, который имеет 27 потоковые мультипроцессоры.

Функция KernelcomputeForces - это медленный удар, занимающий около 95% времени, и мне было интересно, могу ли я что-то сделать дляоптимизировать мой код.

Насколько я вижу, я оптимизировал для памяти пространство-локальность и запись в память.Где-то в документах CUDA говорится, что совместно используемая память быстрее, но я не понимаю, как я могу это использовать.Я разделил работу на блоки 16 с потоками 512 в каждом, но для меня это не совсем понятно.

Пожалуйста, помогите и спасибо за прочтение.

n   is number of bodies

gm  is the gpu mass pointer

gpx is the gpu position x pointer

gpy is the gpu position y pointer

gpz is the gpu position z pointer

gfx is the gpu force x pointer

gfy is the gpu force y pointer

gfz is the gpu force z pointer

Соответствующая функция ядра

__global__ void KernelcomputeForces( unsigned int n, float* gm, float* gpx, float* gpy, float* gpz, float* gfx, float* gfy, float* gfz ){
    int tid = blockDim.x * blockIdx.x + threadIdx.x;
    int numThreads = blockDim.x * gridDim.x;

    float GRAVITY = 0.00001f;

    //compare all with all
    for( unsigned int ia=tid; ia<n; ia+=numThreads ){
        float lfx = 0.0f;
        float lfy = 0.0f;
        float lfz = 0.0f;

        for( unsigned int ib=0; ib<n; ib++ ){
            //compute distance
            float dx = ( gpx[ib] - gpx[ia]);
            float dy = ( gpy[ib] - gpy[ia] );
            float dz = ( gpz[ib] - gpz[ia] );
            //float distance = sqrt( dx*dx + dy*dy + dz*dz );
            float distanceSquared = dx*dx + dy*dy + dz*dz;

            //prevent slingshots and division by zero
            //distance += 0.1f;
            distanceSquared += 0.01f;

            //calculate gravitational magnitude between the bodies
            //float magnitude = GRAVITY * ( gm[ia] * gm[ib] ) / ( distance * distance * distance * distance );
            float magnitude = GRAVITY * ( gm[ia] * gm[ib] ) / ( distanceSquared );

            //calculate forces for the bodies
            //magnitude times direction
            lfx += magnitude * ( dx );
            lfy += magnitude * ( dy );
            lfz += magnitude * ( dz );
        }

        //stores local memory to global memory
        gfx[ia] = lfx;
        gfy[ia] = lfy;
        gfz[ia] = lfz;
    }
}

extern void GPUcomputeForces( unsigned int n, float* gm, float* gpx, float* gpy, float* gpz, float* gfx, float* gfy, float* gfz ){  
    dim3 gridDim( 16, 1, 1 ); //specifys how many blocks in three possible dimensions
    dim3 blockDim( 512, 1, 1 ); //threads per block
    KernelcomputeForces<<<gridDim, blockDim>>>( n, gm, gpx, gpy, gpz, gfx, gfy, gfz );
}

Answer 1

@ talonmies уже ответили на этот вопрос, показывая, как разделяемая память помогает повысить производительность.Таким образом, добавить нечего.Я просто хочу предоставить полный код с различными шагами по оптимизации, включая мозаику , разделяемую память и операции перемешивания , и показать некоторые результаты проведенного тестированияна Kepler K20c.

Приведенный ниже код, сформированный вокруг того, что представлено @talonmies, имеет 5 функций ядра, а именно

KernelcomputeForces: без оптимизации

KernelcomputeForces_Shared: использование плиток в массах источника и общей памяти

KernelcomputeForces_Tiling: использование плиток в массах назначения

KernelcomputeForces_Tiling_Shared: использует мозаику как в исходной, так и в целевой массе и использует общую память

KernelcomputeForces_Tiling_Shuffle: то же, что и выше, но вместо этого используются операции тасованияобщей памяти.

Возможно, по сравнению с тем, что уже опубликовано в литературе ( Быстрое моделирование N-тела с CUDA ) и тем, что уже доступно в виде кодов (см. ответы выше и * 1039).* Марк Харрис 'GitHub N-body page , лАст ядро ​​это единственная новая вещь.Но я немного поиграл с N-body и нашел полезным опубликовать этот ответ, потенциально полезный для следующих пользователей.

Вот код

#include <stdio.h>

#define GRAVITATIONAL_CONST 6.67*1e−11
#define SOFTENING 1e-9f

/********************/
/* CUDA ERROR CHECK */
/********************/
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
    if (code != cudaSuccess) 
    {
        fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
        if (abort) exit(code);
    }
}

/**********/
/* iDivUp */
/**********/
int iDivUp(int a, int b) { return ((a % b) != 0) ? (a / b + 1) : (a / b); }

/*******************/
/* KERNEL FUNCTION */
/*******************/
template<int BLOCKSIZE>
__global__
void KernelcomputeForces(float* m_d, float* x_d, float* y_d, float* z_d, float* fx_d, float* fy_d, float* fz_d, unsigned int N)
{
    int tid = blockDim.x * blockIdx.x + threadIdx.x;

    if (tid < N) {

        float invDist, invDist3;

        // --- Initialize register accumulators for forces
        float fx_temp = 0.0f, fy_temp = 0.0f, fz_temp = 0.0f;

        // --- Move target particle data to registers
        float x_reg = x_d[tid], y_reg = y_d[tid], z_reg = z_d[tid], m_reg = m_d[tid];

        // --- Interact all with all    
        for(unsigned int ib=0; ib<N; ib++) {

            // --- Compute relative distances
            float dx = (x_d[ib] - x_reg);
            float dy = (y_d[ib] - y_reg);
            float dz = (z_d[ib] - z_reg);
            float distanceSquared = dx*dx + dy*dy + dz*dz;

            // --- Prevent slingshots and division by zero
            distanceSquared += SOFTENING;

            float invDist = rsqrtf(distanceSquared);
            float invDist3 = invDist * invDist * invDist;

            // --- Calculate gravitational magnitude between the bodies
            float magnitude = GRAVITATIONAL_CONST * (m_reg * m_d[ib]) * invDist3;

            // --- Calculate forces for the bodies: magnitude times direction
            fx_temp += magnitude*dx;
            fy_temp += magnitude*dy;
            fz_temp += magnitude*dz;
        }

        // --- Stores local memory to global memory
        fx_d[tid] = fx_temp;
        fy_d[tid] = fy_temp;
        fz_d[tid] = fz_temp;
    }
}

/**********************************/
/* KERNEL FUNCTION: SHARED MEMORY */
/**********************************/
template<int BLOCKSIZE>
__global__
void KernelcomputeForces_Shared(float* m_d, float* x_d, float* y_d, float* z_d, float* fx_d, float* fy_d, float* fz_d, unsigned int N)
{
    int tid = blockDim.x * blockIdx.x + threadIdx.x;

    if (tid < N) {

        float invDist, invDist3;

        __shared__ float x_sh[BLOCKSIZE], y_sh[BLOCKSIZE], z_sh[BLOCKSIZE], m_sh[BLOCKSIZE];

        // --- Initialize register accumulators for forces
        float fx_temp = 0.0f, fy_temp = 0.0f, fz_temp = 0.0f;

        // --- Move target particle data to registers
        float x_reg = x_d[tid], y_reg = y_d[tid], z_reg = z_d[tid], m_reg = m_d[tid];

        // --- Interact all with all    
        for(unsigned int ib=0; ib<N; ib+=BLOCKSIZE) {

            // --- Loading data to shared memory
            x_sh[threadIdx.x] = x_d[ib + threadIdx.x];
            y_sh[threadIdx.x] = y_d[ib + threadIdx.x];
            z_sh[threadIdx.x] = z_d[ib + threadIdx.x];
            m_sh[threadIdx.x] = m_d[ib + threadIdx.x];

            __syncthreads();

#pragma unroll
            for(unsigned int ic=0; ic<BLOCKSIZE; ic++) {

                // --- Compute relative distances
                float dx = (x_sh[ic] - x_reg);
                float dy = (y_sh[ic] - y_reg);
                float dz = (z_sh[ic] - z_reg);
                float distanceSquared = dx*dx + dy*dy + dz*dz;

                // --- Prevent slingshots and division by zero
                distanceSquared += SOFTENING;

                float invDist = rsqrtf(distanceSquared);
                float invDist3 = invDist * invDist * invDist;

                // --- Calculate gravitational magnitude between the bodies
                float magnitude = GRAVITATIONAL_CONST * (m_reg * m_sh[ic]) * invDist3;

                // --- Calculate forces for the bodies: magnitude times direction
                fx_temp += magnitude*dx;
                fy_temp += magnitude*dy;
                fz_temp += magnitude*dz;
            }
            __syncthreads();
        }

        // --- Stores local memory to global memory
        fx_d[tid] = fx_temp;
        fy_d[tid] = fy_temp;
        fz_d[tid] = fz_temp;
    }
}

/***************************/
/* KERNEL FUNCTION: TILING */
/***************************/
template<int BLOCKSIZE>
__global__
void KernelcomputeForces_Tiling(float* m_d, float* x_d, float* y_d, float* z_d, float* fx_d, float* fy_d, float* fz_d, unsigned int N)
{
    float invDist, invDist3;

    for (unsigned int tid = blockIdx.x*blockDim.x + threadIdx.x;
                  tid < N;
                  tid += blockDim.x*gridDim.x) {

        // --- Initialize register accumulators for forces
        float fx_temp = 0.0f, fy_temp = 0.0f, fz_temp = 0.0f;

        // --- Move target particle data to registers
        float x_reg = x_d[tid], y_reg = y_d[tid], z_reg = z_d[tid], m_reg = m_d[tid];

        // --- Interact all with all    
        for(unsigned int ib=0; ib<N; ib++) {

            // --- Compute relative distances
            float dx = (x_d[ib] - x_reg);
            float dy = (y_d[ib] - y_reg);
            float dz = (z_d[ib] - z_reg);
            float distanceSquared = dx*dx + dy*dy + dz*dz;

            // --- Prevent slingshots and division by zero
            distanceSquared += SOFTENING;

            float invDist = rsqrtf(distanceSquared);
            float invDist3 = invDist * invDist * invDist;

            // --- Calculate gravitational magnitude between the bodies
            float magnitude = GRAVITATIONAL_CONST * (m_reg * m_d[ib]) * invDist3;

            // --- Calculate forces for the bodies: magnitude times direction
            fx_temp += magnitude*dx;
            fy_temp += magnitude*dy;
            fz_temp += magnitude*dz;
        }
        __syncthreads();

        // --- Stores local memory to global memory
        fx_d[tid] = fx_temp;
        fy_d[tid] = fy_temp;
        fz_d[tid] = fz_temp;
    }
}

/*******************************************/
/* KERNEL FUNCTION: TILING + SHARED MEMORY */
/*******************************************/
template<int BLOCKSIZE>
__global__
void KernelcomputeForces_Tiling_Shared(float* m_d, float* x_d, float* y_d, float* z_d, float* fx_d, float* fy_d, float* fz_d, unsigned int N)
{
    float invDist, invDist3;

    for (unsigned int tid = blockIdx.x*blockDim.x + threadIdx.x;
                  tid < N;
                  tid += blockDim.x*gridDim.x) {

        __shared__ float x_sh[BLOCKSIZE], y_sh[BLOCKSIZE], z_sh[BLOCKSIZE], m_sh[BLOCKSIZE];

        // --- Initialize register accumulators for forces
        float fx_temp = 0.0f, fy_temp = 0.0f, fz_temp = 0.0f;

        // --- Move target particle data to registers
        float x_reg = x_d[tid], y_reg = y_d[tid], z_reg = z_d[tid], m_reg = m_d[tid];

        // --- Interact all with all    
        for(unsigned int ib=0; ib<N; ib+=BLOCKSIZE) {

            // --- Loading data to shared memory
            x_sh[threadIdx.x] = x_d[ib + threadIdx.x];
            y_sh[threadIdx.x] = y_d[ib + threadIdx.x];
            z_sh[threadIdx.x] = z_d[ib + threadIdx.x];
            m_sh[threadIdx.x] = m_d[ib + threadIdx.x];

            __syncthreads();

#pragma unroll
            for(unsigned int ic=0; ic<BLOCKSIZE; ic++) {

                // --- Compute relative distances
                float dx = (x_sh[ic] - x_reg);
                float dy = (y_sh[ic] - y_reg);
                float dz = (z_sh[ic] - z_reg);
                float distanceSquared = dx*dx + dy*dy + dz*dz;

                // --- Prevent slingshots and division by zero
                distanceSquared += SOFTENING;

                float invDist = rsqrtf(distanceSquared);
                float invDist3 = invDist * invDist * invDist;

                // --- Calculate gravitational magnitude between the bodies
                float magnitude = GRAVITATIONAL_CONST * (m_reg * m_sh[ic]) * invDist3;

                // --- Calculate forces for the bodies: magnitude times direction
                fx_temp += magnitude*dx;
                fy_temp += magnitude*dy;
                fz_temp += magnitude*dz;
            }
            __syncthreads();
        }

        // --- Stores local memory to global memory
        fx_d[tid] = fx_temp;
        fy_d[tid] = fy_temp;
        fz_d[tid] = fz_temp;
    }
}

/************************************************/
/* KERNEL FUNCTION: TILING + SHUFFLE OPERATIONS */
/************************************************/
__global__
oid KernelcomputeForces_Tiling_Shuffle(float* m_d, float* x_d, float* y_d, float* z_d, float* fx_d, float* fy_d, float* fz_d, unsigned int N)
{
    float invDist, invDist3;

    const int laneid = threadIdx.x & 31;

    for (unsigned int tid = blockIdx.x*blockDim.x + threadIdx.x;
                  tid < N;
                  tid += blockDim.x*gridDim.x) {

        // --- Initialize register accumulators for forces
        float fx_temp = 0.0f, fy_temp = 0.0f, fz_temp = 0.0f;

        // --- Move target particle data to registers
        float x_reg = x_d[tid], y_reg = y_d[tid], z_reg = z_d[tid], m_reg = m_d[tid];

        // --- Interact all with all    
        for(unsigned int ib=0; ib<N; ib+=32) {

            // --- Loading data to shared memory
            float x_src = x_d[ib + laneid];
            float y_src = y_d[ib + laneid];
            float z_src = z_d[ib + laneid];
            float m_src = m_d[ib + laneid];

#pragma unroll 32
            for(unsigned int ic=0; ic<32; ic++) {

                // --- Compute relative distances
                float dx = (__shfl(x_src, ic) - x_reg);
                float dy = (__shfl(y_src, ic) - y_reg);
                float dz = (__shfl(z_src, ic) - z_reg);
                float distanceSquared = dx*dx + dy*dy + dz*dz;

                // --- Prevent slingshots and division by zero
                distanceSquared += SOFTENING;

                float invDist = rsqrtf(distanceSquared);
                float invDist3 = invDist * invDist * invDist;

                // --- Calculate gravitational magnitude between the bodies
                float magnitude = GRAVITATIONAL_CONST * (m_reg * __shfl(m_src, ic)) * invDist3;

                // --- Calculate forces for the bodies: magnitude times direction
                fx_temp += magnitude*dx;
                fy_temp += magnitude*dy;
                fz_temp += magnitude*dz;
            }
            __syncthreads();
        }

        // --- Stores local memory to global memory
        fx_d[tid] = fx_temp;
        fy_d[tid] = fy_temp;
        fz_d[tid] = fz_temp;
    }
}

/*****************************************/
/* WRAPPER FUNCTION FOR GPU CALCULATIONS */
/*****************************************/
template<int BLOCKSIZE>
float GPUcomputeForces(float* m_d, float* x_d, float* y_d, float* z_d, float* fx_d, float* fy_d, float* fz_d, unsigned int N) {  
    float time;

    dim3 grid(iDivUp(N,BLOCKSIZE), 1, 1);       // --- Specifys how many blocks in three possible dimensions
    dim3 block(BLOCKSIZE, 1, 1);                // --- Threads per block

    cudaEvent_t start, stop;
    gpuErrchk(cudaEventCreate(&start));
    gpuErrchk(cudaEventCreate(&stop));

    gpuErrchk(cudaEventRecord(start, 0));
    KernelcomputeForces_Shared<BLOCKSIZE><<<grid, block>>>(m_d, x_d, y_d, z_d, fx_d, fy_d, fz_d, N);
    //KernelcomputeForces_Tiling<BLOCKSIZE><<<grid, block>>>(m_d, x_d, y_d, z_d, fx_d, fy_d, fz_d, N);
    //KernelcomputeForces<BLOCKSIZE><<<grid, block>>>(m_d, x_d, y_d, z_d, fx_d, fy_d, fz_d, N);
    gpuErrchk(cudaPeekAtLastError());
    gpuErrchk(cudaDeviceSynchronize());

    gpuErrchk(cudaEventRecord(stop, 0));
    gpuErrchk(cudaEventSynchronize(stop));
    gpuErrchk(cudaEventElapsedTime(&time, start, stop));

    return time;
}

/*****************************************/
/* WRAPPER FUNCTION FOR GPU CALCULATIONS */
/*****************************************/
template<int GRIDSIZE, int BLOCKSIZE>
float GPUcomputeForces_Tiling(float* m_d, float* x_d, float* y_d, float* z_d, float* fx_d, float* fy_d, float* fz_d, unsigned int N) {  
    float time;

    dim3 grid(GRIDSIZE, 1, 1);      // --- Specifys how many blocks in three possible dimensions
    dim3 block(BLOCKSIZE, 1, 1);    // --- Threads per block

    cudaEvent_t start, stop;
    gpuErrchk(cudaEventCreate(&start));
    gpuErrchk(cudaEventCreate(&stop));

    gpuErrchk(cudaEventRecord(start, 0));
    //KernelcomputeForces_Tiling<BLOCKSIZE><<<grid, block>>>(m_d, x_d, y_d, z_d, fx_d, fy_d, fz_d, N);
    KernelcomputeForces_Tiling_Shuffle<<<grid, block>>>(m_d, x_d, y_d, z_d, fx_d, fy_d, fz_d, N);
    gpuErrchk(cudaPeekAtLastError());
    gpuErrchk(cudaDeviceSynchronize());

    gpuErrchk(cudaEventRecord(stop, 0));
    gpuErrchk(cudaEventSynchronize(stop));
    gpuErrchk(cudaEventElapsedTime(&time, start, stop));

    return time;
}

/********************/
/* CPU CALCULATIONS */
/********************/
void CPUcomputeForces(float* m_h, float* x_h, float* y_h, float* z_h, float* fx_h, float* fy_h, float* fz_h, unsigned int N) {  

    for (unsigned int i=0; i<N; i++) {

        float invDist, invDist3;

        float fx_temp = 0.0f, fy_temp = 0.0f, fz_temp = 0.0f;

        // --- Interact all with all    
        for(unsigned int ib=0; ib<N; ib++) {

            // --- Compute relative distances
            float dx = (x_h[ib] - x_h[i]);
            float dy = (y_h[ib] - y_h[i]);
            float dz = (z_h[ib] - z_h[i]);
            float distanceSquared = dx*dx + dy*dy + dz*dz;

            // --- Prevent slingshots and division by zero
            distanceSquared += SOFTENING;

            float invDist = 1.f / sqrtf(distanceSquared);
            float invDist3 = invDist * invDist * invDist;

            // --- Calculate gravitational magnitude between the bodies
            float magnitude = GRAVITATIONAL_CONST * (m_h[i] * m_h[ib]) * invDist3;

            // --- Calculate forces for the bodies: magnitude times direction
            fx_temp += magnitude*dx;
            fy_temp += magnitude*dy;
            fz_temp += magnitude*dz;

        }

        // --- Stores local memory to global memory
        fx_h[i] = fx_temp;
        fy_h[i] = fy_temp;
        fz_h[i] = fz_temp;
    }
}

/********/
/* MAIN */
/********/
int main(void)
{
    const int N = 16384;

    size_t gsize = sizeof(float) * size_t(N);

    float * g[10], * _g[7];

    for(int i=0; i<7; i++) gpuErrchk( cudaMalloc((void **)&_g[i], gsize)); 

    for(int i=0; i<10; i++) g[i] = (float *)malloc(gsize);

    for(int i=0; i<N; i++)
        for(int j=0; j<4; j++)
            *(g[j]+i) = (float)rand();

    for(int i=0; i<4; i++) gpuErrchk(cudaMemcpy(_g[i], g[i], gsize, cudaMemcpyHostToDevice)); 

    // --- Warm up to take context establishment time out.
    GPUcomputeForces<512>(_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N);
    //GPUcomputeForces_Tiling<32,512>(_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N);

    // --- Bench runs
    printf("1024: %f\n",    GPUcomputeForces<512>   (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("512:  %f\n",    GPUcomputeForces<512>   (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("256:  %f\n",    GPUcomputeForces<256>   (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("128:  %f\n",    GPUcomputeForces<128>   (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("64:   %f\n",    GPUcomputeForces<64>    (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("32:   %f\n",    GPUcomputeForces<32>    (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));

    printf("16, 1024: %f\n",    GPUcomputeForces_Tiling<16,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("8,  1024: %f\n",    GPUcomputeForces_Tiling<8,512>  (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("4,  1024: %f\n",    GPUcomputeForces_Tiling<4,512>  (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("32, 512: %f\n",     GPUcomputeForces_Tiling<32,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("16, 512: %f\n",     GPUcomputeForces_Tiling<16,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("8,  512: %f\n",     GPUcomputeForces_Tiling<8,512>  (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("64, 256:  %f\n",    GPUcomputeForces_Tiling<64,256> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("32, 256:  %f\n",    GPUcomputeForces_Tiling<32,256> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("16, 256:  %f\n",    GPUcomputeForces_Tiling<16,256> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("128,128:  %f\n",    GPUcomputeForces_Tiling<128,128>(_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("64, 128:  %f\n",    GPUcomputeForces_Tiling<64,128> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("32, 128:  %f\n",    GPUcomputeForces_Tiling<32,128> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("256,64:   %f\n",    GPUcomputeForces_Tiling<256,64> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("128,64:   %f\n",    GPUcomputeForces_Tiling<128,64> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("64, 64:   %f\n",    GPUcomputeForces_Tiling<64,64>  (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("512,32:   %f\n",    GPUcomputeForces_Tiling<512,32> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("256,32:   %f\n",    GPUcomputeForces_Tiling<512,32> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
    printf("128,32:   %f\n",    GPUcomputeForces_Tiling<512,32> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));

    for(int i=8; i<10; i++) gpuErrchk(cudaMemcpy(g[i], _g[i-3], gsize, cudaMemcpyDeviceToHost)); 

    CPUcomputeForces(g[0],g[1],g[2],g[3],g[4],g[5],g[6],N);

    for(int i=0; i<N; i++)
        for(int j=8; j<10; j++) {

            if (abs(*(g[j]+i) - *(g[j-3]+i)) > 0.001f) {
                printf("Error at %i, %i; GPU = %f; CPU = %f\n",i, (j-8), *(g[j-3]+i), *(g[j]+i));
                return;
            }
        }

    printf("Test passed!\n");

    cudaDeviceReset();

    return 0;
}

Вот некоторые результаты дляN = 16384.

KernelcomputeForces: оптимально BLOCKSIZE = 128;t = 10.07ms.

KernelcomputeForces_Shared: оптимально BLOCKSIZE = 128;t = 7.04ms.

KernelcomputeForces_Tiling: оптимально BLOCKSIZE = 128;t = 11.14ms.

KernelcomputeForces_Tiling_Shared: оптимально GRIDSIZE = 64;оптимальный BLOCKSIZE = 256;t = 7.20ms.

KernelcomputeForces_Tiling_Shuffle: оптимально GRIDSIZE = 128;оптимальный BLOCKSIZE = 128;t = 6.84ms.

Кажется, что очень незначительные операции деформирования деформации улучшают производительность по сравнению с общей памятью.

Answer 2

Прежде чем делать что-либо еще, попробуйте запустить больше блоков. Данный блок когда-либо работает только на одном SM - используя только 16 блоков, вы гарантируете, что около 40% емкости графического процессора будет простаивать. Несколько, кратное 27, должно быть оптимальным числом блоков на вашем GTX260-216. Вы также можете обнаружить, что уменьшение количества потоков в блоке не повлияет на производительность, так что вы можете сохранить примерно одинаковый объем работы для каждого потока, но просто сделайте это с достаточным количеством блоков, чтобы охватить весь SM в GPU.

EDIT : Просто чтобы проиллюстрировать это, рассмотрим этот маленький тестовый комплект для вашего ядра:

template<int blocksize, int gridsize>
extern float GPUcomputeForces( unsigned int n, float* gm, float* gpx, float* gpy, float* gpz, float* gfx, float* gfy, float* gfz ){  
    float time;

    dim3 gridDim( gridsize, 1, 1 ); //specifys how many blocks in three possible dimensions
    dim3 blockDim( blocksize, 1, 1 ); //threads per block

    cudaEvent_t start, stop;
    errchk( cudaEventCreate(&start) );
    errchk( cudaEventCreate(&stop) );

    errchk( cudaEventRecord(start, 0) );
    KernelcomputeForces<<<gridDim, blockDim>>>( n, gm, gpx, gpy, gpz, gfx, gfy, gfz );
    rterrchk;

    errchk( cudaEventRecord(stop, 0) );
    errchk( cudaEventSynchronize(stop) );
    errchk( cudaEventElapsedTime(&time, start, stop) );

    return time;
}

int main(void)
{
    const int n = 16384;
    size_t gsize = sizeof(float) * size_t(n);

    float * g[4], * _g[7];

    errchk( cudaSetDevice(1) );  // GTX 275

    for(int i=0; i<7; i++) 
        errchk( cudaMalloc((void **)&_g[i], gsize) ); 

    for(int i=0; i<4; i++)
        g[i] = (float *)malloc(gsize);

    for(int i=0; i<n; i++)
    for(int j=0; j<4; j++)
    *(g[j]+i) = (float)drand48();

    for(int i=0; i<4; i++) 
        errchk( cudaMemcpy(_g[i], g[i], gsize, cudaMemcpyHostToDevice) ); 

    // Warm up to take context establishment time out.
    GPUcomputeForces<16,512>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]);

    // Bench runs
    printf("(1,1)@(512,1,1): %f\n", GPUcomputeForces<1,512>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
    printf("(8,1)@(512,1,1): %f\n", GPUcomputeForces<8,512>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
    printf("(16,1)@(512,1,1): %f\n", GPUcomputeForces<16,512>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
    printf("(30,1)@(256,1,1): %f\n", GPUcomputeForces<30,256>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
    printf("(60,1)@(128,1,1): %f\n", GPUcomputeForces<60,128>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
    printf("(120,1)@(64,1,1): %f\n", GPUcomputeForces<120,64>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
    printf("(240,1)@(32,1,1): %f\n", GPUcomputeForces<240,32>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );

    cudaThreadExit();

    return 0;
}

Когда я запускаю это, я получаю следующие времена выполнения на стандартном GTX 275 (все времена в миллисекундах):

(1,1)@(512,1,1): 1087.107910
(8,1)@(512,1,1): 135.582458
(16,1)@(512,1,1): 67.876671
(30,1)@(256,1,1): 54.881279
(60,1)@(128,1,1): 35.261280
(120,1)@(64,1,1): 36.316288
(240,1)@(32,1,1): 39.870495

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

Answer 3

Совместно используемая память будет полезной оптимизацией в ядре такого типа - она ​​позволяет объединять чтения положений и масс частиц, что на GT200 будет очень важно.Я обнаружил, что это примерно в два раза быстрее, чем ваша версия (запущенная с вашими 16384 частицами с использованием 128 блоков из 128 потоков):

template<int blocksize>
__global__
void KernelcomputeForces1( unsigned int n1, float* gm, float* gpx,
                float* gpy, float* gpz, float* gfx, float* gfy, float* gfz )
{
    __shared__ float lgpx[blocksize], lgpy[blocksize],
               lgpz[blocksize], lgm[blocksize];

    const float GRAVITY = 0.00001f;

    //compare all with all
    float lfx = 0.0f, lfy = 0.0f, lfz = 0.0f;
    int ia = blockDim.x * blockIdx.x + threadIdx.x;
    float lgpx0 = gpx[ia], lgpy0 = gpy[ia],
          lgpz0 = gpz[ia], lgm0 = gm[ia];

    for( unsigned int ib=0; ib<n1; ib+=blocksize ){

        lgpx[threadIdx.x] = gpx[ib + threadIdx.x];
        lgpy[threadIdx.x] = gpy[ib + threadIdx.x];
        lgpz[threadIdx.x] = gpz[ib + threadIdx.x];
        lgm[threadIdx.x] = gm[ib + threadIdx.x];
        __syncthreads();

#pragma unroll
        for(unsigned int ic=0; ic<blocksize; ic++) {

            //compute distance
            float dx = ( lgpx[ic] - lgpx0 );
            float dy = ( lgpy[ic] - lgpy0 );
            float dz = ( lgpz[ic] - lgpz0 );
            float distanceSquared = dx*dx + dy*dy + dz*dz;

            //prevent slingshots and division by zero
            distanceSquared += 0.01f;

            //calculate gravitational magnitude between the bodies
            float magnitude = GRAVITY * ( lgm0 * lgm[ic] )
                    / ( distanceSquared );

            //calculate forces for the bodies
            //magnitude times direction
            lfx += magnitude * ( dx );
            lfy += magnitude * ( dy );
            lfz += magnitude * ( dz );
        }
        __syncthreads();
    }

    //stores local memory to global memory
    gfx[ia] = lfx;
    gfy[ia] = lfy;
    gfz[ia] = lfz;
}

Вам нужно будет сделать что-то немного другое для числа частицкоторые выходят за пределы кратного размера блока, вероятно, вторая строфа, которая не будет развернута.Следите за возможностью расхождения деформации с помощью вызовов __syncthreads (), которые могут привести к зависанию ядра, если вы не будете осторожны.

Answer 4

Профилировщик cuda расскажет вам информацию о занятости, условном ветвлении, глобальном использовании памяти (пропуски кэша или несообщенные чтения, в зависимости от версии cuda) и многое другое. Это важная часть оптимизации вашего ядра.

Answer 5

Я хочу предоставить отдельный, дополнительный ответ на ответ, приведенный выше, который, я думаю, будет полезен для следующих пользователей.

Способ дальнейшей оптимизации проблемы N-тела - прибегнуть к дереву.подходы, основанные на Barnes-Hut, которые были распараллелены в

M. Burtscher, K. Pingali, "An Efficient CUDA Implementation of the Tree-Based Barnes Hut n-Body

«Алгоритм», GPU Computing Gems, Emerald Edition.

и чья реализация CUDA доступна для загрузки по адресу LonestarGPU .

Сравнивая вышеупомянутые ядра на основе разделяемых и случайных чисел с упомянутой реализацией на Kepler K20c, я получил следующие результаты (времена в мс)

N           Barnes-Hut          Shared          Shuffle
16384       19.1                5.7             7.2
32768       46.9                25.5            21.7
65536       107.7               102.6           74.6
131072      255.1               355             296.2
262144      548.3               1408.8          1108.9
524288      1246                5434            4688
1048576     2674                21544           18632
2097152     5664                85980           74454

Обратите внимание, что такой анализ недоступен в оригинальной публикации, указанной выше.