Matrix inverse usng linear system solver через cublas, cublasCreate exception или что-то еще

Question

Я пытаюсь инвертировать матрицу с помощью решения линейного уравнения через библиотеку CUDA cublas.

Исходное уравнение в виде:

Ax = B = I

I - единичная матрица А - матрица Я пытаюсь обратного х - (надеюсь) обратное (А) матрица

Я хотел бы выполнить разложение LU, получают следующее:

LUx = в

L - это нижний треугольник матрицы U - это верхняя треугольная матрица

Вот хороший пример, который может показать, что я пытаюсь сделать:

http://www.youtube.com/watch?v=dza5JTvMpzk

Ради обсуждения , давайте предположим, что у меня уже есть разложение LU A. (A = LU) Я пытаюсь найти обратную последовательно в два этапа:

Ax = B = I

LUx = B = I

1-й шаг: Объявляет у:

Ly = I

решения 1-ый линейное уравнение

2-й шаг, решить следующее линейное уравнение

Ux = y

найти X = обратный (A) - Bingo! @!

---

Albeit, сейчас я понятия не имею, что я делаю неправильно здесь ...

Там может быть два предположения, либо я не использую cublas правильно или cublas бросает исключение для нет причина.продолжайте читать :)

---

Смотрите мой код прилагается, Eсть код MATLAB в конце:

---

 #include "cuda_runtime.h" 
    #include "device_launch_parameters.h" 

    #include <stdio.h> 


    //#include "cublas.h" 
    #include "cublas_v2.h" 


int main() 
{ 

cudaError_t cudaStat; 
cublasStatus_t stat; 
cublasHandle_t handle; 

//cublasInit(); 

stat = cublasCreate(&handle); 
if (stat != CUBLAS_STATUS_SUCCESS) 
{ 
    printf ("CUBLAS initialization failed\n"); 
    return -1; 
} 

unsigned int size = 3; 
// Allowcate device matrixes 
float* p_h_LowerTriangle; 
float* p_d_LowerTriangle; 

p_h_LowerTriangle = new float[size*size]; 
p_h_LowerTriangle[0] = 1.f; 
p_h_LowerTriangle[1] = 0.f; 
p_h_LowerTriangle[2] = 0.f; 
p_h_LowerTriangle[3] = 2.56f; 
p_h_LowerTriangle[4] = 1.f; 
p_h_LowerTriangle[5] = 0.f; 
p_h_LowerTriangle[6] = 5.76f; 
p_h_LowerTriangle[7] = 3.5f; 
p_h_LowerTriangle[8] = 1.f; 

float* p_h_UpperTriangle; 
float* p_d_UpperTriangle; 

p_h_UpperTriangle = new float[size*size]; 
p_h_UpperTriangle[0] = 25.f; 
p_h_UpperTriangle[1] = 5.f; 
p_h_UpperTriangle[2] = 1.f; 
p_h_UpperTriangle[3] = 0.f; 
p_h_UpperTriangle[4] = -4.8f; 
p_h_UpperTriangle[5] = -1.56f; 
p_h_UpperTriangle[6] = 0.f; 
p_h_UpperTriangle[7] = 0.f; 
p_h_UpperTriangle[8] = 0.7f; 

float* p_h_IdentityMatrix; 
float* p_d_IdentityMatrix; 

p_h_IdentityMatrix = new float[size*size]; 
p_h_IdentityMatrix[0] = 1.f; 
p_h_IdentityMatrix[1] = 0.f; 
p_h_IdentityMatrix[2] = 0.f; 
p_h_IdentityMatrix[3] = 0.f; 
p_h_IdentityMatrix[4] = 1.f; 
p_h_IdentityMatrix[5] = 0.f; 
p_h_IdentityMatrix[6] = 0.f; 
p_h_IdentityMatrix[7] = 0.f; 
p_h_IdentityMatrix[8] = 1.f; 

cudaMalloc ((void**)&p_d_LowerTriangle, size*size*sizeof(float)); 
cudaMalloc ((void**)&p_d_UpperTriangle, size*size*sizeof(float)); 
cudaMalloc ((void**)&p_d_IdentityMatrix, size*size*sizeof(float)); 


float* p_d_tempMatrix; 
cudaMalloc ((void**)&p_d_tempMatrix, size*size*sizeof(float)); 


stat = cublasSetMatrix(size,size,sizeof(float),p_h_LowerTriangle,size,p_d_LowerTriangle,size); 
stat = cublasSetMatrix(size,size,sizeof(float),p_h_UpperTriangle,size,p_d_UpperTriangle,size); 
stat = cublasSetMatrix(size,size,sizeof(float),p_h_IdentityMatrix,size,p_d_IdentityMatrix,size); 

cudaDeviceSynchronize(); 

// 1st phase - solve Ly = b 
const float alpha = 1.f; 

// Function solves the triangulatr linear system with multiple right hand sides, function overrides b as a result 
stat = cublasStrsm(handle, 
    cublasSideMode_t::CUBLAS_SIDE_LEFT, 
    cublasFillMode_t::CUBLAS_FILL_MODE_LOWER, 
    cublasOperation_t::CUBLAS_OP_N, 
    CUBLAS_DIAG_UNIT, 
    size, 
    size, 
    &alpha, 
    p_d_LowerTriangle, 
    size, 
    p_d_IdentityMatrix, 
    size); 

//////////////////////////////////// 
// TODO: printf 1st phase the results 
cudaDeviceSynchronize(); 

cudaMemcpy(p_h_IdentityMatrix,p_d_IdentityMatrix,size*size*sizeof(float),cudaMemcpyDeviceToHost); 
cudaMemcpy(p_h_UpperTriangle,p_d_UpperTriangle,size*size*sizeof(float),cudaMemcpyDeviceToHost); 
cudaMemcpy(p_h_LowerTriangle,p_d_LowerTriangle,size*size*sizeof(float),cudaMemcpyDeviceToHost); 

stat =cublasGetMatrix(size,size,sizeof(float),p_d_IdentityMatrix,size,p_h_IdentityMatrix,size); 
stat =cublasGetMatrix(size,size,sizeof(float),p_d_UpperTriangle,size,p_h_UpperTriangle,size); 
stat =cublasGetMatrix(size,size,sizeof(float),p_d_LowerTriangle,size,p_h_LowerTriangle,size); 

//////////////////////////////////// 

// 2nd phase - solve Ux = b 
stat = cublasStrsm(handle, 
    cublasSideMode_t::CUBLAS_SIDE_LEFT, 
    cublasFillMode_t::CUBLAS_FILL_MODE_UPPER, 
    cublasOperation_t::CUBLAS_OP_N, 
    CUBLAS_DIAG_NON_UNIT, 
    size, 
    size, 
    &alpha, 
    p_d_UpperTriangle, 
    size, 
    p_d_IdentityMatrix, 
    size); 

// TODO print the results 
cudaDeviceSynchronize(); 

//////////////////////////////////// 
cudaMemcpy(p_h_IdentityMatrix,p_d_IdentityMatrix,size*size*sizeof(float),cudaMemcpyDeviceToHost); 
cudaMemcpy(p_h_UpperTriangle,p_d_UpperTriangle,size*size*sizeof(float),cudaMemcpyDeviceToHost); 
cudaMemcpy(p_h_LowerTriangle,p_d_LowerTriangle,size*size*sizeof(float),cudaMemcpyDeviceToHost); 
//////////////////////////////////// 

cublasDestroy(handle); 
//cublasShutdown(); 

cudaFree(p_d_LowerTriangle); 
cudaFree(p_d_UpperTriangle); 
cudaFree(p_d_IdentityMatrix); 

free(p_h_LowerTriangle); 
free(p_h_UpperTriangle); 
free(p_h_IdentityMatrix); 


return 0; 
} 

---

Matlab код - работает идеально

---

/////////////////////////// MATLAB

clear all 

    UpperMatrix = [ 25 5 1 ; 0 -4.8 -1.56 ; 0 0 0.7 ] 

    LowerMatrix = [1 0 0 ; 2.56 1 0 ; 5.76 3.5 1 ] 

    IdentityMatrix = eye(3) 

    % 1st phase solution 
    otps1.LT = true; 
    y = linsolve(LowerMatrix,IdentityMatrix,otps1) 

    %2nd phase solution 
    otps2.UT = true; 
    y = linsolve(UpperMatrix,y,otps2) 

//////////////////////////////// MATLAB выход

UpperMatrix =

25,0000 5,0000 1,0000 0 -4,8000 -1,5600 0 0 0,7000

LowerMatrix =

1.0000   0   0 
2.5600 1.0000   0 
5.7600 3.5000 1.0000 

IdentityMatrix =

1  0  0 
0  1  0 
0  0  1 

у =

1.0000   0   0 

-2,5600 1,0000 0 3,2000 -3,5000 1,0000

у =

0.0476 -0.0833 0.0357 

-0,9524 1,4167 -0,4643 4,5714 -5,0000 1,4286

UpperMatrix =

25,0000 5,0000 1,0000 0 -4,8000 -1,5600 0 0 0,7000

LowerMatrix =

1.0000   0   0 
2.5600 1.0000   0 
5.7600 3.5000 1.0000 

IdentityMatrix =

1  0  0 
0  1  0 
0  0  1 

у =

1.0000   0   0 

-2,5600 1,0000 0 3,2000 -3,5000 1.0000

Inverse_LU_UT

UpperMatrix =

25,0000 5,0000 1,0000 0 -4,8000 -1,5600 0 0 0,7000

LowerMatrix =

1.0000   0   0 
2.5600 1.0000   0 
5.7600 3.5000 1.0000 

IdentityMatrix =

1  0  0 
0  1  0 
0  0  1 

у =

1.0000   0   0 

-2,5600 1,0000 0 3,2000 -3,5000 1,0000

у =

0.0476 -0.0833 0.0357 

-0,9524 1,4167 -0,4643 4,5714 -5,0000 1,4286

---

Я понятия не имею, что я делаю неправильно здесь, я подозреваю, что операция cublasCreate,

Каждый раз, когда я попал в команду:

cublasStatus_t stat; 
    cublasHandle_t handle; 
    stat = cublasCreate(&handle); 

стат и ручкой переменные действительны.

НО выход VS10 несколько сообщений об ошибках укажите следующее:

первый шанс исключение в 0x ... Microsoft C++ исключение cudaError_enum в ячейке памяти 0x ...

Некоторые могут сказать, что это внутренняя ошибка cublas сообщение, обрабатываемое самой библиотекой.

Я Ted грустно сейчас,

Помощь будет очень признателен :)

Answer 1

Вы знаете, что, cuBlas хранит матрицы в столбцах способа, но Matlab и использование C-строка из основных способов.

Измените порядок инициализации и запуска. Результат должен быть хорошим.