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

Вы упираетесь в пределы точности float. На самом деле не следует ожидать, что он будет иметь точность более 5 знаков после запятой, и некоторые шаблоны вычислений могут привести к тому, что он будет иметь меньшую точность. На самом деле, результат CUBLAS численно ближе к правильно округленному результату, чем ваша реализация ЦП. Это легко доказать. Все, что нам нужно сделать, это выполнить операцию суммирования вашего хоста, используя double, и мы увидим, что получили другой результат:

$ cat t1784.cu
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cublas_v2.h>
#define HANDLE_ERROR(x) x
void creatematrix(float *out, int nx, int ny)
{
    float ctr = 0;
    for (int i = 0; i < nx; ++i) {
        for (int j = 0; j < ny; ++j) {
            out[j * nx + i] = ctr/1E+5;
            ctr = ctr + 1;
        }
    }
}

float add_arr_val(float *im, int N)
{
    float tmp = 0;
    for (int i = 0; i < N; ++i)
        tmp += im[i];

    float out = tmp;
    return out;
}

double add_arr_val_dbl(float *im, int N)
{
    double tmp = 0;
    for (int i = 0; i < N; ++i)
        tmp += (double)(im[i]);

    return tmp;
}


__global__ void init_ones(float *d_in, int N)
{
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    if (i < N)
    {
        d_in[i] = 1.0;
    }
}

int main()
{
    // Define matrix size (using flattened array for most operations)
    int nx = 512;       // row size
    int ny = 512;       // column size
    int N = nx * ny;    // total size of flattened array

    // CPU section ========================================
    float *M; M = (float*)malloc(N * sizeof(float));    // create array pointer and allocate memory
    creatematrix(M, nx, ny);                            // create a test matrix of size nx * ny
    float cpu_out = add_arr_val(M, N);                  // CPU function
    double cpu_dbl_out = add_arr_val_dbl(M, N);
    // GPU and cuBLAS section ==============================
    float *d_M;
    HANDLE_ERROR(cudaMalloc(&d_M, N * sizeof(float)));
    HANDLE_ERROR(cudaMemcpy(d_M, M, N * sizeof(float), cudaMemcpyHostToDevice));    // copy original array M to device as d_M

    // create array of all ones, size N for dot product
    float *d_ones;
    cudaMalloc(&d_ones, N * sizeof(float));

    // Max potential blocksize
    int minGridSize, blockSize, gridSize;
    cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, init_ones, 0, N);
    gridSize = (N + blockSize - 1) / blockSize;
    init_ones << <gridSize, blockSize >> > (d_ones, N);     // kernel launch to generate array of all 1's
    cudaDeviceSynchronize();

    float blas_out;                                                         // output on host variable
    cublasHandle_t handle;  cublasCreate(&handle);                          // initialize CUBLAS context
    cublasSdot(handle, N, d_M, 1, d_ones, 1, &blas_out);                    // Perform cublas single-precision dot product of (d_M . d_ones)
    cudaDeviceSynchronize();

    //print output from cpu and gpu sections
    printf("native output = %f\n", cpu_out);
    printf("native double output = %f\n", cpu_dbl_out);
    printf("cublas output = %f\n", blas_out);

    cublasDestroy(handle);
    free(M);
    cudaFree(d_M);
    cudaFree(d_ones);
}
$ nvcc -o t1784 t1784.cu -lcublas
$ ./t1784
native output        = 343590.437500
native double output = 343596.072960
cublas output        = 343596.062500
$

Причина, по которой вывод cublas на самом деле ближе (*), заключается в том, что он выполняет добавления в другом порядке, чем ваш код хоста float. Он работает в блоках потоков и суммирует частичные суммы перед созданием окончательного результата.

В качестве примечания: нет необходимости использовать l со спецификатором формата %f printf. Он уже предназначен для правильной передачи форматов double и float.

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

(*) Тем не менее, вы не должны предполагать, что это всегда так, хотя я считаю, что метод частичных сумм, как правило, более надежен, чем метод чистой текущей суммы (это просто личное мнение, а не вопрос, который я доказал или что хочу поспорить) Однако в любом случае ошибка зависит от данных. Мы можем построить конкретную последовательность данных, которая сделает метод скользящей суммы очень хорошим с точки зрения точности. Для выполнения большого суммирования, где важен наивысший уровень точности, вам, вероятно, следует использовать максимально возможную точность. Кроме того, вы можете прочитать о суммировании Кахана в статье, ссылка на которую приведена выше.