#cuda #simulation
#cuda #Симуляция
Вопрос:
Я разрабатываю алгоритм с N-телом в 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 );
}
Комментарии:
1. Я согласен, я не вижу немедленного способа использовать разделяемую память в этом коде, поскольку каждый поток обращается ко всей глобальной памяти.
Ответ №1:
@talonmies уже ответил на этот вопрос, показывая, как разделяемая память полезна для повышения производительности. Итак, добавить больше нечего. Я просто хочу предоставить полный код с различными этапами оптимизации, включая разбиение на листы, разделяемую память и операции перетасовки, и показать некоторые результаты тестирования, выполненного на Kepler K20c.
Приведенный ниже код, созданный на основе того, что представлено @talonmies, имеет 5 функции ядра, а именно
KernelcomputeForces : оптимизация отсутствует
KernelcomputeForces_Shared : использует тайлинг в массивах исходных текстов и совместно используемой памяти
KernelcomputeForces_Tiling : использует тайлинг в целевых массивах
KernelcomputeForces_Tiling_Shared : использует как тайлинги в массивах источника и назначения, так и общую память.
KernelcomputeForces_Tiling_Shuffle : то же, что и выше, но использует операции случайного перемещения вместо разделяемой памяти.
Возможно, по сравнению с тем, что уже опубликовано в литературе (быстрое моделирование N-Body с помощью CUDA) и тем, что уже доступно в виде кодов (см. Приведенные Выше ответы и N-body страницу Марка Харриса на GitHub, последнее ядро является единственной новой вещью. Но я немного поиграл с 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 %dn", 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 amp; 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(amp;start));
gpuErrchk(cudaEventCreate(amp;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(amp;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(amp;start));
gpuErrchk(cudaEventCreate(amp;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(amp;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 **)amp;_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: %fn", GPUcomputeForces<512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("512: %fn", GPUcomputeForces<512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("256: %fn", GPUcomputeForces<256> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("128: %fn", GPUcomputeForces<128> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("64: %fn", GPUcomputeForces<64> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("32: %fn", GPUcomputeForces<32> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("16, 1024: %fn", GPUcomputeForces_Tiling<16,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("8, 1024: %fn", GPUcomputeForces_Tiling<8,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("4, 1024: %fn", GPUcomputeForces_Tiling<4,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("32, 512: %fn", GPUcomputeForces_Tiling<32,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("16, 512: %fn", GPUcomputeForces_Tiling<16,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("8, 512: %fn", GPUcomputeForces_Tiling<8,512> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("64, 256: %fn", GPUcomputeForces_Tiling<64,256> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("32, 256: %fn", GPUcomputeForces_Tiling<32,256> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("16, 256: %fn", GPUcomputeForces_Tiling<16,256> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("128,128: %fn", GPUcomputeForces_Tiling<128,128>(_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("64, 128: %fn", GPUcomputeForces_Tiling<64,128> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("32, 128: %fn", GPUcomputeForces_Tiling<32,128> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("256,64: %fn", GPUcomputeForces_Tiling<256,64> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("128,64: %fn", GPUcomputeForces_Tiling<128,64> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("64, 64: %fn", GPUcomputeForces_Tiling<64,64> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("512,32: %fn", GPUcomputeForces_Tiling<512,32> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("256,32: %fn", GPUcomputeForces_Tiling<512,32> (_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6],N));
printf("128,32: %fn", 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 = %fn",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 ; optimal BLOCKSIZE = 256 ; t = 7.20ms .
KernelcomputeForces_Tiling_Shuffle : оптимальный GRIDSIZE = 128 ; optimal BLOCKSIZE = 128 ; t = 6.84ms .
Похоже, что операции искажения в случайном порядке очень незначительно улучшают производительность по сравнению с общей памятью.
Ответ №2:
Общая память будет полезной оптимизацией в ядре такого типа — она позволяет объединять считывания положений и масс частиц, что на 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(), которые могут привести к зависанию ядра, если вы не будете осторожны.
Комментарии:
1. интересно, хотя, если у меня 27 мультипроцессоров, не было бы оптимальным использовать 27 потоков?
2. Нет. Каждому мультипроцессору на графическом процессоре класса gt200 требуется не менее 192 потоков только для покрытия задержки конвейера команд. Вы хотите, чтобы количество блоков было кратно 27 — любому блоку, выполняемому только на одном мультипроцессоре, при этом каждый блок имеет поток, кратный 32 потокам, что является базовой детализацией планирования графического процессора.
3. хорошо, спасибо всем, я получил все, что хотел знать сейчас. далее, моделирование частиц воды: D
Ответ №3:
Прежде чем делать что-либо еще, попробуйте запустить больше блоков. Данный блок всегда выполняется только на одном SM — используя только 16 блоков, вы гарантируете, что около 40% мощности графического процессора будет простаивать. Некоторое кратное 27 должно быть оптимальным количеством блоков на вашем GTX260-216. Вы также можете обнаружить, что уменьшение количества потоков на блок не повредит производительности, так что вы можете выполнять примерно одинаковый объем работы на поток, но просто делайте это с достаточным количеством блоков, чтобы покрыть все SM в графическом процессоре.
РЕДАКТИРОВАТЬ: Просто чтобы проиллюстрировать суть, рассмотрим этот небольшой тестовый пакет для вашего ядра:
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(amp;start) );
errchk( cudaEventCreate(amp;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(amp;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 **)amp;_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): %fn", GPUcomputeForces<1,512>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
printf("(8,1)@(512,1,1): %fn", GPUcomputeForces<8,512>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
printf("(16,1)@(512,1,1): %fn", GPUcomputeForces<16,512>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
printf("(30,1)@(256,1,1): %fn", GPUcomputeForces<30,256>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
printf("(60,1)@(128,1,1): %fn", GPUcomputeForces<60,128>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
printf("(120,1)@(64,1,1): %fn", GPUcomputeForces<120,64>(n,_g[0],_g[1],_g[2],_g[3],_g[4],_g[5],_g[6]) );
printf("(240,1)@(32,1,1): %fn", 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
Ie: Запуск по крайней мере такого количества блоков, сколько MP на карте, имеет решающее значение для повышения производительности, даже при использовании блоков очень малого размера.
Комментарии:
1. то, что вы только что попросили меня сделать, было просто неправильно. изменение с 16 на> 27 уменьшило мой fps до 16. изменение с 512 на> 256 уменьшило мой fps до 10
2. @ColacX — вам захочется попробовать МНОЖЕСТВО конфигураций ядра. Различные конфигурации оказывают существенное влияние и обычно имеют положительные стороны.
3. @ColacX — Я очень удивлен этим. Мой собственный сравнительный анализ вашего ядра на аналогичной плате (GTX 275) ясно показывает, что время выполнения ядра с по крайней мере таким количеством блоков, как у вашей карты SM, сокращает время выполнения ядра почти вдвое. Именно так, как я и предсказывал. Смотрите мою правку для получения дополнительной информации.
4.
code (16,1)@(512,1,1): 36.050186 ms (32,1)@(256,1,1): 35.151672 ms (48,1)@(256,1,1): 34.227173 ms (64,1)@(128,1,1): 31.856787 ms (27,1)@(512,1,1): 40.098667 ms (54,1)@(64,1,1): 60.681877 ms (54,1)@(128,1,1): 40.320881 ms (54,1)@(256,1,1): 36.719494 ms (54,1)@(512,1,1): 39.685223 ms (81,1)@(64,1,1): 50.386776 ms (81,1)@(128,1,1): 31.069546 ms (81,1)@(256,1,1): 36.522266 ms (81,1)@(512,1,1): 39.731533 ms (108,1)@(64,1,1): 40.844597 ms (108,1)@(128,1,1): 35.667061 ms (108,1)@(256,1,1): 35.841301 ms (108,1)@(512,1,1): 39.554295 msкратные 27, похоже, не так сильно улучшают5. Вы действительно уверены, что используете этот код на GTX 260-216? Эти цифры на стендах выглядят совершенно иначе, чем на GTX 275, который я использовал. Кроме того, GTX 260 должен быть примерно на 15-20% медленнее, чем 275 (более низкий уровень шейдеров и тактовая частота памяти), но ваши результаты показывают примерно на 50% более высокую производительность при наименьшем размере блока. У меня большой опыт настройки кода для этих карт — что-то здесь вообще не имеет смысла.
Ответ №4:
Я хочу предоставить отдельный дополнительный ответ к приведенному выше, который, я думаю, будет полезен следующим пользователям.
Способ дальнейшей оптимизации задачи из N элементов — прибегнуть к древовидным подходам, подобным подходу Барнса-Хата, который был распараллелен в
M. Burtscher, K. Pingali, "An Efficient CUDA Implementation of the Tree-Based Barnes Hut n-Body
Алгоритм», GPU Computing Gems, издание Emerald.
и чья реализация 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
Пожалуйста, обратите внимание, что такой анализ недоступен в оригинальной публикации, указанной выше.
Ответ №5:
Профилировщик cuda предоставит вам информацию о занятости, условном ветвлении, глобальном использовании памяти (пропуски кэша или несогласованные чтения, в зависимости от версии cuda) и многом другом. Это важная часть оптимизации вашего ядра.
Комментарии:
1. о, не знал, что есть один, я проверю это, хотя, держу пари, в профилировщике я не вижу ничего очевидного. все еще могли бы пригодиться некоторые советы по изменению кода