Обратная матрица cuBLAS намного медленнее, чем MATLAB

В моем текущем проекте я пытаюсь вычислить обратную большую (n > 2000) матрицу с помощью cuBLAS. Выполняется обратный расчет, но по какой-то причине время расчета значительно меньше, чем при выполнении в MATLAB.

Я приложил образец вычисления, выполненного на случайных матрицах с использованием моей реализации на любом языке, а также результаты производительности.

Мы будем очень признательны за любую помощь или предложения о том, что может быть причиной этого замедления.

Заранее спасибо.

Сравнение

cuBLAS против MATLAB

N = 500: cuBLAS ~ 0,130 с, MATLAB ~ 0,066 с -> ~ в 1,97 раза медленнее

N = 1000: cuBLAS ~ 0,898 с, MATLAB ~ 0,311 с -> ~ в 2,89 раза медленнее

N = 2000: cuBLAS ~ 6,667 с, MATLAB ~ 0,659 с -> ~ в 10,12 раза медленнее

N = 4000: cuBLAS ~ 51,860 с, MATLAB ~ 4,296 с -> ~ в 12,07 раз медленнее

Код C++

#include <string>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <conio.h>

#define CUDA_CALL(res, str) { if (res != cudaSuccess) { printf("CUDA Error : %s : %s %d : ERR %s\n", str, __FILE__, __LINE__, cudaGetErrorName(res)); } }
#define CUBLAS_CALL(res, str) { if (res != CUBLAS_STATUS_SUCCESS) { printf("CUBLAS Error : %s : %s %d : ERR %d\n", str, __FILE__, __LINE__, int(res)); } }

static cudaEvent_t cu_TimerStart;
static cudaEvent_t cu_TimerStop;

void d_CUDATimerStart(void)
{
    CUDA_CALL(cudaEventCreate(&cu_TimerStart), "Failed to create start event!");
    CUDA_CALL(cudaEventCreate(&cu_TimerStop), "Failed to create stop event!");

    CUDA_CALL(cudaEventRecord(cu_TimerStart), "Failed to record start event!");
}

float d_CUDATimerStop(void)
{
    CUDA_CALL(cudaEventRecord(cu_TimerStop), "Failed to record stop event!");

    CUDA_CALL(cudaEventSynchronize(cu_TimerStop), "Failed to synch stop event!");

    float ms;

    CUDA_CALL(cudaEventElapsedTime(&ms, cu_TimerStart, cu_TimerStop), "Failed to elapse events!");

    CUDA_CALL(cudaEventDestroy(cu_TimerStart), "Failed to destroy start event!");
    CUDA_CALL(cudaEventDestroy(cu_TimerStop), "Failed to destroy stop event!");

    return ms;
}

float* d_GetInv(float* L, int n)
{
    cublasHandle_t cu_cublasHandle;
    CUBLAS_CALL(cublasCreate(&cu_cublasHandle), "Failed to initialize cuBLAS!");

    float** adL;
    float** adC;
    float* dL;
    float* dC;
    int* dLUPivots;
    int* dLUInfo;

    size_t szA = n * n * sizeof(float);

    CUDA_CALL(cudaMalloc(&adL, sizeof(float*)), "Failed to allocate adL!");
    CUDA_CALL(cudaMalloc(&adC, sizeof(float*)), "Failed to allocate adC!");
    CUDA_CALL(cudaMalloc(&dL, szA), "Failed to allocate dL!");
    CUDA_CALL(cudaMalloc(&dC, szA), "Failed to allocate dC!");
    CUDA_CALL(cudaMalloc(&dLUPivots, n * sizeof(int)), "Failed to allocate dLUPivots!");
    CUDA_CALL(cudaMalloc(&dLUInfo, sizeof(int)), "Failed to allocate dLUInfo!");

    CUDA_CALL(cudaMemcpy(dL, L, szA, cudaMemcpyHostToDevice), "Failed to copy to dL!");
    CUDA_CALL(cudaMemcpy(adL, &dL, sizeof(float*), cudaMemcpyHostToDevice), "Failed to copy to adL!");
    CUDA_CALL(cudaMemcpy(adC, &dC, sizeof(float*), cudaMemcpyHostToDevice), "Failed to copy to adC!");

    d_CUDATimerStart();

    CUBLAS_CALL(cublasSgetrfBatched(cu_cublasHandle, n, adL, n, dLUPivots, dLUInfo, 1), "Failed to perform LU decomp operation!");
    CUDA_CALL(cudaDeviceSynchronize(), "Failed to synchronize after kernel call!");

    CUBLAS_CALL(cublasSgetriBatched(cu_cublasHandle, n, (const float **)adL, n, dLUPivots, adC, n, dLUInfo, 1), "Failed to perform Inverse operation!");
    CUDA_CALL(cudaDeviceSynchronize(), "Failed to synchronize after kernel call!");

    float timed = d_CUDATimerStop();

    printf("cublas inverse in: %.5f ms.\n", timed);

    float* res = (float*)malloc(szA);

    CUDA_CALL(cudaMemcpy(res, dC, szA, cudaMemcpyDeviceToHost), "Failed to copy to res!");

    CUDA_CALL(cudaFree(adL), "Failed to free adL!");
    CUDA_CALL(cudaFree(adC), "Failed to free adC!");
    CUDA_CALL(cudaFree(dL), "Failed to free dL!");
    CUDA_CALL(cudaFree(dC), "Failed to free dC!");
    CUDA_CALL(cudaFree(dLUPivots), "Failed to free dLUPivots!");
    CUDA_CALL(cudaFree(dLUInfo), "Failed to free dLUInfo!");

    CUBLAS_CALL(cublasDestroy(cu_cublasHandle), "Failed to destroy cuBLAS!");

    return res;
}

int main()
{
    int n = 1000;
    float* L = (float*)malloc(n * n * sizeof(float));
    for(int i = 0; i < n * n; i++)
        L[i] = ((float)rand()/(float)(RAND_MAX));

    float* inv = d_GetInv(L, n);

    printf("done.");
    _getch();

    return 0;
}

Код MATLAB

A = rand(1000);
tic
X = inv(A);
toc

Системная информация:

Видеокарта: GTX 780 3 ГБ

ЦП: i7-4790S @ 3,20 ГГц


person at1012    schedule 09.06.2016    source источник
comment
Используемый вами метод не предназначен для быстрой обратной матрицы для больших матриц в CUBLAS. Вы читали документацию? Указывает: эта функция предназначен для использования с матрицами небольших размеров, где существенным фактором являются накладные расходы на запуск.   -  person Robert Crovella    schedule 09.06.2016
comment
@RobertCrovella, приношу свои извинения, я, должно быть, проглядел эту строчку. Есть ли какие-либо другие средства для вычисления большой матричной инверсии, которые стоит изучить с помощью библиотек cuda? или реализация Matlab более быстрый путь?   -  person at1012    schedule 09.06.2016


Ответы (1)


Как сказал @RobertCrovella, вы не должны использовать пакетные API-интерфейсы для небольших матриц для одной большой матричной инверсии.

В основном вы можете использовать тот же метод, что и в вашем коде, но с непакетной версией getrf() и getri() для максимальной производительности для большой матрицы.

Для getrf() вы можете найти его здесь.

http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-getrf

Для getri(), хотя инструментарий CUDA не предоставляет getri() для решения AX=I, где A является LU-фактором getrf(), он предоставляет getrs() для решения AX=B. Все, что вам нужно сделать, это установить B=I перед вызовом getrs().

http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-getrs

person kangshiyin    schedule 09.06.2016
comment
При реализации вашего решения я генерирую идентификационную матрицу B в памяти устройства, что является узким местом в скорости и использовании памяти. Я действительно недавно задавал вопрос об этом. Знаете ли вы, можно ли его сгенерировать быстрее, чем наивный цикл for? Или если эту проблему можно решить с меньшим использованием памяти? - person avgn; 17.06.2018
comment
@avgn Я не думаю, что настройка B=I должна быть узким местом в скорости. Вы можете указать более подробную информацию в своем вопросе. - person kangshiyin; 25.08.2018