MATLAB производит другой результат, чем CUBLAS + Kernel

Question

У меня есть следующий MATLAB код:

[N, d] = size(X); % data size and dimensions 

R = rand(d,dt); % Form a random matrix with elements in [0,1] 

% Random projection 
Y = X * R; 

w=720; % hashing step 

b = w * rand(dt,1); 

% Compute the hash codes of the data 
binId = floor(bsxfun(@plus, Y, b')/w); 

, и я попытался сделать его параллельно с использованием CUBLAS и ядра следующим образом:

__global__ void compute(const int N,const int dt,const int w,const float *old, int *newt){ 
    int col = blockDim.y * blockIdx.y + threadIdx.y; 
    int row = blockDim.x * blockIdx.x + threadIdx.x; 
    int id = row+N*col; 
    if(row<N && col<dt){ 
     newt[id]=(floor)(old[id]/w); 
    } 
} 

void gpu_blas_mmul(cublasHandle_t handle, const float *A, const float *B, float *C, const int m, const int k, const int n, const float bet) { 
    int lda=m,ldb=k,ldc=m; 
    const float alf = 1.0; 
    const float *alpha = &alf; 
    const float *beta = &bet; 

    // Do the actual multiplication and addition 
    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc); 
} 

float *d_R, *d_RX, *d_B_row; 
int *d_H; 

thrust::device_vector<float> d_X(h_X, h_X + N * d); 

cudaMalloc(&d_R,d * dt * sizeof(float)); 
cudaMemcpy(d_R,h_R,d * dt * sizeof(float),cudaMemcpyHostToDevice); 

cudaMalloc(&d_B_row,dt * sizeof(float)); 
cudaMemcpy(d_B_row,h_B_row,dt * sizeof(float),cudaMemcpyHostToDevice); 

cudaMalloc(&d_RX,N * dt * sizeof(float)); 
cudaMalloc(&d_H,N * dt * sizeof(int)); 

//-------------------------CuBLAS----------------------- 

cublasHandle_t handle; 
cublasCreate(&handle); 

thrust::device_vector<float> d_B_col(N, w); 

gpu_blas_mmul(handle, thrust::raw_pointer_cast(&d_B_col[0]), d_B_row, d_RX, N, 1, dt,0.0); 

gpu_blas_mmul(handle, thrust::raw_pointer_cast(&d_X[0]), d_R, d_RX, N, d, dt, 1.0); 

cublasDestroy(handle); 

//-----------------------Kernel---------------------------- 
dim3 blockSize(BLOCK_SIZE, BLOCK_SIZE,1); 
int linGrid1 = (int)ceil(N/(float)BLOCK_SIZE); 
int linGrid2 = (int)ceil(dt/(float)BLOCK_SIZE); 
dim3 gridSize(linGrid1,linGrid2,1); 
compute<<<gridSize, blockSize>>>(N, dt, w, d_RX, d_H); 

В h_X, h_R и h_B_row Я сохранил (в основном порядке) X, R и b, созданные MATLAB. Используемый мною набор данных - ANN_SIFT1M от http://corpus-texmex.irisa.fr/

Для примерно 10000 значений полученные результаты являются точно такими же, но когда я пытаюсь использовать 50000 значений, например, есть некоторые различия, которые становятся все больше и больше по мере увеличения количества значений.

Любое представление о том, что я делаю неправильно?

Answer 1

Ваш код MATLAB использует точность с двумя точками, поэтому результат является более точным. В отличие от этого, ядро ​​CUDA, которое вы предоставили, использует точную точность, тип float и, следовательно, дает менее точный результат. И как обычно, когда речь идет о проблемах с одиночной и двойной точками, проблема только ухудшается, когда вы начинаете увеличивать размер ваших входных данных.

Решение должно использовать тип double вместо float.