Числовая ошибка в простом ядре cuda / cublas с использованием определенного ввода
Я работаю с cuda и cublas, и я пытался реализовать простые операции, такие как матричное поэлементное умножение / деление. Я использую только поплавок для моих экспериментов. Я знаю, что наиболее очевидный способ сделать это — написать ядро, подобное этому:
__global__ void mul_elementwise(const unsigned int n, float* source, float* dest, const float value)
{
const unsigned int offset = blockIdx.x * blockDim.x + threadIdx.x;
const unsigned int stride = blockDim.x * gridDim.x;
for (unsigned int i = offset; i < n; i += stride)
{
dest[i] = source[i] * value;
}
}
Это ядро может работать как для умножения, так и для деления (просто используя 1 / x в качестве значения). Но этого также можно достичь, используя библиотеку cublas: предположим, что у нас есть матрица A mxn, хранящаяся в стиле старшей колонки, и скаляр x, а затем в качестве вектора m * n 1s устанавливаются alpha = x или alpha = 1 / x и d_ones, мы можем вызвать и получить тот же результат
cublasSaxpy(cublas_handle, m * n, &alpha, d_ones, 1, A_dev, 1);
Оба метода работают просто отлично, но я сталкиваюсь с несколькими проблемами с определенной матрицей, для которой оба метода не работают. Я выделил эту большую матрицу и построил MCVE доступны Вот (вы можете скомпилировать его с nvcc mcve.cu -lcublas, Как вы можете видеть, результаты в обоих случаях совершенно неверны: результаты хоста совершенно разные, я пытаюсь выяснить, что происходит. Я не вижу никаких ошибок в коде, но, возможно, я должен попытаться использовать double вместо float и посмотреть, что произойдет.
Есть мнения по поводу этой ситуации? Заранее спасибо!
РЕДАКТИРОВАТЬ # 1 Я попытался использовать double, но ничего не изменилось, если я использую cublasDaxpy, в то же время он отлично работает с собственным ядром. Я думаю, что значения слишком малы, поэтому одной точности с плавающей точкой недостаточно.
Решение
Интересный MCVE. Разве не было бы возможно уменьшить ваш вектор до нескольких элементов? Разве нельзя показать расхождение в расчете только на 1 элемент вектора?
Во всяком случае я вижу несколько проблем.
-
Ваше ядро реализует следующую функцию:
y=alpha*x, Но SAXPY инвентарьy=alpha*x+y, Сейчас еслиyначинался как (все) ноль, тогда эти два были бы одинаковыми. Но это не то, что у вас есть:CUBLAS Your Kernel --------------------------- alpha: alpha alpha x: 1 ahost (ahost is your huge data array) y: ahost -Итак, ваше ядро вычисляет
y=alpha * ahost, но ваш вызов CUBLAS вычисляетy = alpha*1 + ahost, Я не ожидал бы такого же результата от них, в общем. -
Ваш анализ ошибки кажется ошибочным в нескольких отношениях. Во-первых, вы вычисляете абсолютную ошибку в
floatпеременная (число, которое всегда будет положительным, поскольку это абсолютное значение), но затем вы сравниваете его с отрицательный число:float diff = abs(host[i]-dev[i]); ... if (diff > (-1e12))не так ли
ifпроверка всегда будет правдой? Возможно, вы имели в виду1e-12хотя это все равно будет ошибочным. Поиск фиксированного порога ошибки при сравнении с плавающей запятой должен быть масштабирован до размера сравниваемых чисел.floatколичества содержат только около 6-7 точных десятичных цифр. (И суммирование этих ошибок также хлопотно.)
Вот полный код, который исправил вышеуказанные проблемы и выдает ошибку с нулевой суммой для всех сравнений (хост<-> ядро и хост<-> cublas):
static float array[] = {0x00000000,
0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0xB58DA1CF,0xB50D2FEC,0x34A48536,0xB4A1D5BC,0x358E1345,0x35943AAC,0xB5983F40,0xB43628BB,0xB4A95348,0xB4DB751C,0xB50C8D1A,0xB3EFCBB5,0x3552B8CD,0x3538A167,0x358FDE0D,0xB4D54CE9,0xB5D29BB7,0xB4A234EE,0x346EF2F4,0x35B5D9F2,0xB40F1487,0x3554BC20,0x33FD9466,0xB536D37D,0xB3C2E594,0xB59DA581,0x3584FC87,0x34438F09,0x35D293CB,0xB4FBB002,0xB59F41E9};
#include <iostream>
#include <stdio.h>
#include <cublas_v2.h>
#include <assert.h>
#define TOL 0.0001
typedef unsigned int u32;
#define GET_STRIDE() u32(blockDim.x * gridDim.x)
#define GET_OFFSET() u32(blockIdx.x * blockDim.x + threadIdx.x)
inline
cudaError_t checkCuda(cudaError_t result)
{
#if defined(DEBUG) || defined(_DEBUG)
if (result != cudaSuccess) {
fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
assert(result == cudaSuccess);
}
#endif
return result;
}
__global__ void div_elementwise(const u32 n, float* source, float* dest, const float value)
{
for (u32 i = GET_OFFSET(); i < n; i += GET_STRIDE())
{
dest[i] = source[i] * value;
}
}
float check_eq(float* dev, float* host, u32 len)
{
float sum = 0.0f;
for (u32 i = 0; i < len; ++i)
{
if (dev[i]!=host[i])
{
//printf("diff %d %f %f\n", i, dev[i], host[i]);
//break;
float diff = abs((host[i]-dev[i])/host[i]);
sum += diff;
if (diff > (TOL))
printf("diff %d %f\n", i, diff);
}
}
printf("%f\n", sum);
return sum;
}
void div_host(float* a, float v, u32 len)
{
for (u32 i = 0; i < len; ++i)
{
a[i]=a[i]*v;
}
}
int main()
{
u32 len = sizeof(array)/sizeof(float);
printf("array len = %d\n", len);
for (int i =0; i < len; i++) if (isnan(array[i])) {printf("nan value at %d\n",i); return -1;}
float* adev, *adevcublas, *d_zero;
float* ahost = (float*) malloc(len * sizeof(float));
checkCuda(cudaMalloc(&adev, len * sizeof(float)));
checkCuda(cudaMalloc(&adevcublas, len * sizeof(float)));
checkCuda(cudaMalloc(&d_zero, len * sizeof(float)));
memcpy(ahost, &array[0], len * sizeof(float));
checkCuda(cudaMemcpy(adev, ahost, len * sizeof(float), cudaMemcpyHostToDevice));
checkCuda(cudaMemcpy(adevcublas, ahost, len * sizeof(float), cudaMemcpyHostToDevice));
checkCuda(cudaMemset(d_zero, 0, len*sizeof(float)));
float alpha = 1/2494.f;
printf("%f\n", alpha);
div_host(ahost, alpha, len);
u32 tb = 256;
div_elementwise<<<((len + tb - 1) / tb),tb>>>(len, adev, adev, alpha);
float* r = (float*) malloc(len * sizeof(float));
checkCuda(cudaMemcpy(r, adev, len * sizeof(float), cudaMemcpyDeviceToHost));
check_eq(r,ahost,len);cublasHandle_t ch;
cublasCreate(&ch);
float* r0 = (float*) malloc(len * sizeof(float));
cublasStatus_t stat = cublasSaxpy(ch, len, &alpha, adevcublas, 1, d_zero, 1);
if (stat != CUBLAS_STATUS_SUCCESS) {std::cout << "CUBLAS error: " << (int)stat << std::endl; return 1;}
checkCuda(cudaMemcpy(r0, d_zero, len * sizeof(float), cudaMemcpyDeviceToHost));check_eq(r0,ahost,len);
free(r);
free(r0);
free(ahost);
cudaFree(adev);
return 0;
}
Другие решения