Образец cuSolver с cuSolverDnDgetrf не работает

ХОРОШО. Мои руки пачкаются кодом, который я взял из примеров cuSolver.
У меня мало опыта работы с C ++, но мне как-то удалось взять то, что мне нужно, из исходного кода.

Проблема в том, когда я пытаюсь это выполнить; в соответствии с рекомендациями справочного руководства, которое я компилирую:

nvcc -c att3_cus_lu.cpp -I/usr/local/cuda-8.0/targets/x86_64-linux/include
g++ -fopenmp -o res.out att3_cus_lu.o -L/usr/local/cuda/lib64 -lcudart -lcublas -lcusolver -lcusparse

Нет проблем до здесь; результат, который я получаю, всегда один и тот же:

step 1: set matrix (A) and right hand side vector (b)
step 2: prepare data on device
step 3: solve A*x = b
Error: LU factorization failed
INFO_VALUE = 2
timing: LU =   0.000208 sec
step 4: evaluate residual
|b - A*x| = NAN
|A| = 1.000000E+00
|x| = NAN
|b - A*x|/(|A|*|x|) = NAN

Какую бы матрицу или вектор b я ни поставил, результат будет одинаковым; код не может разложить матрицу
Я показал INFO_VALUE, значение которого всегда равно 2 при каждом выполнении; это ценность info переменная, запрошенная для функции cuSolverDnDgetrf (). В справочном руководстве cuSolver сказано:

Эта функция вычисляет факторизацию LU матрицы P ×A = LU
где A — матрица m × n, P — матрица перестановок, L — нижняя треугольная матрица с диагональю единиц, а U — верхняя треугольная матрица.
Если факторизация LU не удалась, то есть матрица A (U) является единственной, выходной параметр devInfo = i указывает U (i, i) = 0.

В приведенном ниже коде я поместил одинаковую матрицу, поэтому никакие особые матрицы не работают.

Вот весь код; шаблон для кодов cuda main () повторяется: задайте переменные хоста, cudaMemcpy их на устройство, выполните их на устройстве с помощью функций cuda, cudaMemcpy их верните на хост, затем продолжайте работу с последовательным кодом, пока не повторитесь.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <assert.h>

#include <cuda_runtime.h>

#include "cublas_v2.h"#include "cusolverDn.h"#include "helper_cuda.h"
#include "helper_cusolver.h"
#ifndef MAX
#define MAX(a,b) (a > b) ? a : b
#endif

void linearSolverLU(
cusolverDnHandle_t handle,
int n,
const double *Acopy,
int lda,
const double *b,
double *x)
{
int bufferSize = 0;
int *info = NULL;
double *buffer = NULL;
double *A = NULL;
int *ipiv = NULL; // pivoting sequence
int h_info = 0;
double start, stop;
double time_solve;

checkCudaErrors(cusolverDnDgetrf_bufferSize(handle, n, n,(double*)Acopy, lda, &bufferSize));

checkCudaErrors(cudaMalloc(&info, sizeof(int)));
checkCudaErrors(cudaMalloc(&buffer, sizeof(double)*bufferSize));
checkCudaErrors(cudaMalloc(&A, sizeof(double)*lda*n));
checkCudaErrors(cudaMalloc(&ipiv, sizeof(int)*n));// prepare a copy of A because getrf will overwrite A with L
checkCudaErrors(cudaMemcpy(A, Acopy, sizeof(double)*lda*n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemset(info, 0, sizeof(int)));

start = second();

checkCudaErrors(cusolverDnDgetrf(handle, n, n, A, lda, buffer, ipiv, info));
checkCudaErrors(cudaMemcpy(&h_info, info, sizeof(int), cudaMemcpyDeviceToHost));

if ( 0 != h_info ){
fprintf(stderr, "Error: LU factorization failed\n");
printf("INFO_VALUE = %d\n",h_info);
}

checkCudaErrors(cudaMemcpy(x, b, sizeof(double)*n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cusolverDnDgetrs(handle, CUBLAS_OP_N, n, 1, A, lda, ipiv, x, n, info));
checkCudaErrors(cudaDeviceSynchronize());
stop = second();

time_solve = stop - start;
fprintf (stdout, "timing: LU = %10.6f sec\n", time_solve);

if (info  ) { checkCudaErrors(cudaFree(info  )); }
if (buffer) { checkCudaErrors(cudaFree(buffer)); }
if (A     ) { checkCudaErrors(cudaFree(A)); }
if (ipiv  ) { checkCudaErrors(cudaFree(ipiv));}

}

//int main (int argc, char *argv[]) !!!
int main(void)
{
cusolverDnHandle_t handle = NULL;
cublasHandle_t cublasHandle = NULL; // used in residual evaluation
cudaStream_t stream = NULL;

int rowsA = 3; // number of rows of A
int colsA = 3; // number of columns of A
int lda = MAX(colsA, rowsA); // leading dimension in dense matrix

double *h_A = NULL; // dense matrix
double *h_x = NULL; // a copy of d_x
double *h_b = NULL; // b = ones(m,1)
double *h_r = NULL; // r = b - A*x, a copy of d_r

double *d_A = NULL; // a copy of h_A
double *d_x = NULL; // x = A \ b
double *d_b = NULL; // a copy of h_b
double *d_r = NULL; // r = b - A*x

// the constants are used in residual evaluation, r = b - A*x
const double minus_one = -1.0;
const double one = 1.0;

double x_inf = 0.0;
double r_inf = 0.0;
double A_inf = 0.0;
int errors = 0;

int i, j, col, row, k;

h_A = (double*)malloc(sizeof(double)*lda*colsA);
h_x = (double*)malloc(sizeof(double)*colsA);
h_b = (double*)malloc(sizeof(double)*rowsA);
h_r = (double*)malloc(sizeof(double)*rowsA);
assert(NULL != h_A);
assert(NULL != h_x);
assert(NULL != h_b);
assert(NULL != h_r);

memset(h_A, 0., sizeof(double)*lda*colsA);

printf("step 1: set matrix (A) and right hand side vector (b)\n");

double mat[9] = {1.,0.,0.,0.,1.,0.,0.,0.,1.};
double bb[3] = {1., 1., 1.}; //RANDOM MATRICES 4 TESTING

for( row =0; row < rowsA ; row++ )
{

for( col = 0; col < colsA ; col++ )
{
h_A[row*lda + col] = mat[row];
}
}

for( k = 0; k < rowsA; k++ ){
h_b[k] = bb[k];
}

checkCudaErrors(cusolverDnCreate(&handle));
checkCudaErrors(cublasCreate(&cublasHandle));
checkCudaErrors(cudaStreamCreate(&stream));

checkCudaErrors(cusolverDnSetStream(handle, stream));
checkCudaErrors(cublasSetStream(cublasHandle, stream));checkCudaErrors(cudaMalloc((void **)&d_A, sizeof(double)*lda*colsA));
checkCudaErrors(cudaMalloc((void **)&d_x, sizeof(double)*colsA));
checkCudaErrors(cudaMalloc((void **)&d_b, sizeof(double)*rowsA));
checkCudaErrors(cudaMalloc((void **)&d_r, sizeof(double)*rowsA));

printf("step 2: prepare data on device\n");
checkCudaErrors(cudaMemcpy(d_A, h_A, sizeof(double)*lda*colsA, cudaMemcpyHostToDevice));
checkCudaErrors(cudaMemcpy(d_b, h_b, sizeof(double)*rowsA, cudaMemcpyHostToDevice));

printf("step 3: solve A*x = b \n");

linearSolverLU(handle, rowsA, d_A, lda, d_b, d_x);

printf("step 4: evaluate residual\n");
checkCudaErrors(cudaMemcpy(d_r, d_b, sizeof(double)*rowsA, cudaMemcpyDeviceToDevice));

// r = b - A*x
checkCudaErrors(cublasDgemm_v2(
cublasHandle,
CUBLAS_OP_N,
CUBLAS_OP_N,
rowsA,
1,
colsA,
&minus_one,
d_A,
lda,
d_x,
rowsA,
&one,
d_r,
rowsA));

checkCudaErrors(cudaMemcpy(h_x, d_x, sizeof(double)*colsA, cudaMemcpyDeviceToHost));
checkCudaErrors(cudaMemcpy(h_r, d_r, sizeof(double)*rowsA, cudaMemcpyDeviceToHost));

x_inf = vec_norminf(colsA, h_x);
r_inf = vec_norminf(rowsA, h_r);
A_inf = mat_norminf(rowsA, colsA, h_A, lda);

printf("|b - A*x| = %E \n", r_inf);
printf("|A| = %E \n", A_inf);
printf("|x| = %E \n", x_inf);
printf("|b - A*x|/(|A|*|x|) = %E \n", r_inf/(A_inf * x_inf));

if (handle) { checkCudaErrors(cusolverDnDestroy(handle)); }
if (cublasHandle) { checkCudaErrors(cublasDestroy(cublasHandle)); }
if (stream) { checkCudaErrors(cudaStreamDestroy(stream)); }

if (h_A) { free(h_A); }
if (h_x) { free(h_x); }
if (h_b) { free(h_b); }
if (h_r) { free(h_r); }

if (d_A) { checkCudaErrors(cudaFree(d_A)); }
if (d_x) { checkCudaErrors(cudaFree(d_x)); }
if (d_b) { checkCudaErrors(cudaFree(d_b)); }
if (d_r) { checkCudaErrors(cudaFree(d_r)); }

cudaDeviceReset();

return 0;
}

Вот и все. Я знаю, что это много вещей; любая помощь будет оценена.
Я действительно не понимаю, что я делаю не так с этими матрицами.

Дайте мне знать, если я должен добавить другую важную информацию.

Обратите внимание: в исходном коде была запрошена внешняя разреженная матрица с расширением .mtx, а затем в функции linearSolverLU она была преобразована в плотную матрицу. Я удалил этот материал, и теперь я прямо пытаюсь решить линейную систему с плотной матрицей.