низкая производительность OpenCL на HD4850

Здравствуйте.

Мне для диплома нужны расчёты на GPU, но вот сделать их действительно быстрыми не получается.

У меня задача моделирования частиц которые взаимодействуют каждая с каждой, т.е. сложность N^2 где N количество частиц.
Демка nVidia nBodyCS на моей видюхе считает 64k точек с fps ~4
Моя же программа считает 20k точек с fps ~1
nBodyCS написана на DirectCompute и взаимодействи между частицами у меня сложнее поэтому напрямую сравнивать конечно нельзя но всё же.

Как сделать быстрее?

вот код:

  1. инициализация буферов
  2. clmemFV = clCreateBuffer(GPUContext, CL_MEM_READ_ONLY, sizeof(cl_float4)*m_paramColculate.countSteps*m_countParticleAddedAtStap, NULL, &err);
  3. clmemFC = clCreateBuffer(GPUContext, CL_MEM_READ_ONLY, sizeof(cl_float4)*m_paramColculate.countSteps*m_countParticleAddedAtStap, NULL, &err);
  4. clmemFA = clCreateBuffer(GPUContext, CL_MEM_READ_ONLY, sizeof(cl_float4)*m_paramColculate.countSteps*m_countParticleAddedAtStap, NULL, &err);
  5. clmemSV = clCreateBuffer(GPUContext, CL_MEM_WRITE_ONLY, sizeof(cl_float4)*m_paramColculate.countSteps*m_countParticleAddedAtStap, NULL, &err);
  6. clmemSC = clCreateBuffer(GPUContext, CL_MEM_WRITE_ONLY, sizeof(cl_float4)*m_paramColculate.countSteps*m_countParticleAddedAtStap, NULL, &err);
  7. clmemSA = clCreateBuffer(GPUContext, CL_MEM_WRITE_ONLY, sizeof                 (cl_float4)*m_paramColculate.countSteps*m_countParticleAddedAtStap, NULL, &err);
  8.  
  9. расчёт:
  10. err += clEnqueueWriteBuffer ( GPUCommandQueue, clmemFV, CL_TRUE, 0, sizeof(cl_float4)*countPart, m_Velocity[m1s], 0, NULL, NULL);
  11. err += clEnqueueWriteBuffer ( GPUCommandQueue, clmemFC, CL_TRUE, 0, sizeof(cl_float4)*countPart, m_Coordinate[m1s], 0, NULL, NULL);
  12. err += clEnqueueWriteBuffer ( GPUCommandQueue, clmemFA, CL_TRUE, 0, sizeof(cl_float4)*countPart, m_Acceleration_m1[m1s], 0, NULL, NULL);
  13.  
  14. err += clSetKernelArg(clTest, 0, sizeof(cl_mem), (void*)&clmemFV);
  15. err += clSetKernelArg(clTest, 1, sizeof(cl_mem), (void*)&clmemFC);
  16. err += clSetKernelArg(clTest, 2, sizeof(cl_mem), (void*)&clmemFA);
  17. err += clSetKernelArg(clTest, 3, sizeof(cl_mem), (void*)&clmemSV);
  18. err += clSetKernelArg(clTest, 4, sizeof(cl_mem), (void*)&clmemSC);
  19. err += clSetKernelArg(clTest, 5, sizeof(cl_mem), (void*)&clmemSA);
  20. err += clSetKernelArg(clTest, 6, sizeof(cl_int), (void*)&countPart);
  21. err += clSetKernelArg(clTest, 7, sizeof(cl_float), (void*)&t);
  22. err += clSetKernelArg(clTest, 8, sizeof(cl_float4), (void*)&VECTOR4(    m_paramColculate.directionE.x*m_paramColculate.E,
  23.                         m_paramColculate.directionE.y*m_paramColculate.E,
  24.                         m_paramColculate.directionE.z*m_paramColculate.E,
  25.                         0.0));
  26. err += clSetKernelArg(clTest, 9, sizeof(cl_float4), (void*)&VECTOR4(    m_paramColculate.directionB.x*m_paramColculate.B,
  27.                         m_paramColculate.directionB.y*m_paramColculate.B,
  28.                         m_paramColculate.directionB.z*m_paramColculate.B,
  29.                         0.0));
  30. err += clSetKernelArg(clTest, 10, sizeof(cl_float), (void*)&k);
  31.  
  32.                         //***********************EXECUTING KERNEL**********************
  33. WorkSize[0] = m_countParticle;
  34. cl_event  events[1];
  35. err += clEnqueueNDRangeKernel(GPUCommandQueue, clTest, 1, NULL, WorkSize, NULL, 0, NULL, &events[0]);
  36. err += clWaitForEvents(1, &events[0]);
  37. err += clReleaseEvent(events[0]);
  38.  
  39.                         //***********************READING RESULTS**********************
  40. err += clEnqueueReadBuffer(GPUCommandQueue, clmemSV, CL_TRUE, 0, sizeof(cl_float4)*countPart, m_Velocity[cs], 0, NULL, NULL);
  41. err += clEnqueueReadBuffer(GPUCommandQueue, clmemSC, CL_TRUE, 0, sizeof(cl_float4)*countPart, m_Coordinate[cs], 0, NULL, NULL);
  42. err += clEnqueueReadBuffer(GPUCommandQueue, clmemSA, CL_TRUE, 0, sizeof(cl_float4)*countPart, m_Acceleration_m1[cs], 0, NULL, NULL);
  43.  
  44. и ядро:
  45. inline float4 F( float4 currentV, float4 E, float4 B)
  46. {
  47.         float4 a;
  48.         float EV = (E.x*currentV.x + E.y*currentV.y + E.z*currentV.z)/8.9875517e+16f;
  49.         float rel = -1.0f*sqrt(1.0f - (currentV.x*currentV.x + currentV.y*currentV.y + currentV.z*currentV.z)/8.9875517e+16f)*1.756e+11f;
  50.  
  51.         a.x = ( E.x + currentV.y*B.z - currentV.z*B.y - currentV.x*EV ) * rel;
  52.         a.y = ( E.y + currentV.z*B.x - currentV.x*B.z - currentV.y*EV ) * rel;
  53.         a.z = ( E.z + currentV.x*B.y - currentV.y*B.x - currentV.z*EV ) * rel;
  54.         return a;
  55. }
  56. inline float4 mul(float4 A, float4 B)
  57. {
  58.         float4 tmp;
  59.         tmp.x = A.y*B.z - A.z*B.y;
  60.         tmp.y = A.z*B.x - A.x*B.z;
  61.         tmp.z = A.x*B.y - A.y*B.x;
  62.         tmp.w = 0.0f;
  63.         return tmp;
  64. }
  65. __kernel void computeRungeKutta2(       __global float4* FV,
  66.                                                                         __global float4* FC,
  67.                                                                         __global float4* FA,
  68.                                                                         __global float4* SV,
  69.                                                                         __global float4* SC,
  70.                                                                         __global float4* SA,
  71.                                                                         const int size,
  72.                                                                         const float t,
  73.                                                                         const float4 E,
  74.                                                                         const float4 B,
  75.                                                                         const float k
  76.                                                                         )
  77. {
  78.  
  79.         int index = get_global_id(0);
  80.         if(index >= size) return;
  81.  
  82.         const float c = 2.99792458e+8f;
  83.         const float PI = 3.141592653589f;
  84.         float4 f0, f1, tmpV, Ecs = {0.0f, 0.0f, 0.0f, 0.0f}, Bcs = {0.0f, 0.0f, 0.0f, 0.0f};
  85.        
  86.         for(int j = 0; j < size; j++)
  87.         {
  88.                 if(index != j)
  89.                 {                      
  90.                         float4 Rj;
  91.                         Rj = FC[index] - FC[j];
  92.  
  93.                         float Rjl = sqrt(Rj.x*Rj.x + Rj.y*Rj.y + Rj.z*Rj.z);
  94.                        
  95.                         float tmp = 4.0f*PI*8.8541878f*pow((Rjl - (Rj.x*FV[j].x + Rj.y*FV[j].y + Rj.z*FV[j].z)/c),3.0f)/(-1.602e-7f*k);
  96.                         float tmp2 = 1.0f - (FV[j].x*FV[j].x + FV[j].y*FV[j].y + FV[j].z*FV[j].z)/(c*c);
  97.                         float4 tmp1 = mul(Rj, mul((Rj - FV[j] * Rjl/c), FA[j]));
  98.                        
  99.                         float4 Ecsj;
  100.                                
  101.                         Ecsj = ((Rj - FV[j] * Rjl/c) * tmp2 + tmp1/(c*c))/ tmp;
  102.                                                                                
  103.                         Ecs += Ecsj;
  104.                         Bcs += mul(Rj,Ecsj)/(c*Rjl);
  105.                 }
  106.         }
  107.         Bcs.w = 0.0f;
  108.         Ecs.w = 0.0f;
  109.                                                
  110.         f0 = F( FV[index], E + Ecs, B + Bcs);
  111.         f1 = F( FV[index] + f0*t, E + Ecs, B + Bcs);
  112.         SA[index] = f0;
  113.  
  114.         tmpV = FV[index]+ (f0 + f1)*t/2.0f;
  115.         SC[index] = FC[index] + tmpV*t;
  116.         SV[index] = tmpV;
  117.  
  118. }

всем спасибо.

Forums: 

Вместо pow() и sqrt()

Вместо pow() и sqrt() использовать native_pow() и native_sqrt(). Вместо умножений и сложений можно использовать mad(). Некоторые операции можно заменить на dot(). Если видеокарта - ATI, лучше все векторизовать. Функцию mul(float4 A, float4 B) можно заменить/перпеписать как A.yzx*B.zxy - A.zxy*B.yzx. Дальше - смотреть профайлером что просаживает производительность.

Ещё, не надо определять PI - уже есть M_PI_F.

Все операции производятся над 3-х компонентными векторами, но используются 4-х компонентные float4. Если это из соображений быстрого чтения из памяти - проверьте, возможно на платформе и так sizeof(cl_float4) == sizeof(float3), и не имеет смысла использовать float4.