cuSolverサンプルがcuSolverDnDgetrfで動作しない

Question

OK。私は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 

私は結果を入れても同じです。コードは行列を因数分解できません。 私はINFO_VALUEを表示しましたが、その値は実行ごとに常に2です。それはcuSolverDnDgetrf（）関数のために要求された変数infoの値です。 cuSolver参照マニュアルにおいては言う：

この関数は、Aは、×n行列AMのn行列P A = L U ×AMのLU分解を計算し、Pは順列行列であり、Lは単位対角の下三角行列、Uは上三角行列です。 LU分解が失敗した場合、すなわち、行列A（U）出力パラメータのDevInfo、単数形である=私はかかわらず、同一のマトリックスを置かれているの下U（I、i）はコードで0

を=示し特異行列は動いていない。

ここにコード全体があります。 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; 
} 

これだけです。私はこれがたくさんあることを知っています。どんな助けもありがとう。 私は実際に私がそれらのマトリックスで間違っているものを得ることはありません。

その他の関連情報を追加する必要がある場合はお知らせください。

NB元のコードでは、拡張子が.mtxの外部スパース行列が要求され、linearSolverLU関数では密行列に変換されました。私はそのことを取り除きましたが、今では代わりに高密度マトリクスを使って線形システムを解くことを直接試みています。

Answer 1

cuSolverに提出して分解する行列が無効です。ライブラリは、行列エントリのエラーを正しく報告しています。

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

h_Aが特異であり、そして、あなたの意図がなかったものを（私は仮定）{ 1, 1, 1, 0, 0, 0, 0, 0, 0 }を含む生成されます。犯人は、このコードです。