расчет собственного значения и собственного вектора с помощью cusolver из cuda 7.0 RC

Я пытаюсь вычислить наибольшую пару собственное значение/собственный вектор с помощью cuSolver, выпущенного в CUDA 7.0 RC. Проблема в том, что я получаю CUSOLVER_INTERNAL_ERROR и не знаю, что с этим делать.

Это моя удобная вещь, используемая для вызова функций cuda/cusparse/cusolver.

// my handy stuff
#define CUDA_CALL(value) do {                                                     \
    cudaError_t _m_cudaStat = value;                                                  \
    if (_m_cudaStat != cudaSuccess) {                                               \
        fprintf(stderr, "Error %s at line %d in file %s\N",                 \
                cudaGetErrorString(_m_cudaStat), __LINE__, __FILE__);       \
        exit(-1);                                                                             \
    } 
} while(0)

#define CUSPARSE_CALL(value) do {                                                                 \
    cusparseStatus_t _m_status = value;                                                             \
    if (_m_status != CUSPARSE_STATUS_SUCCESS){                                                      \
        fprintf(stderr, "Error %d at line %d in file %s\N", (int)_m_status, __LINE__, __FILE__);    \
    exit(-5);                                                                                     \
    }                                                                                             \
} while(0)

#define CUSOLVER_CALL(value) do {                                                                 \
    cusolverStatus_t _m_status = value;                                                             \
    if (_m_status != CUSOLVER_STATUS_SUCCESS){                                                      \
        fprintf(stderr, "Error %d at line %d in file %s\N", (int)_m_status, __LINE__, __FILE__);    \
    exit(-5);                                                                                     \
    }                                                                                             \
} while(0)

И это мой код

#include "cusparse.h"
#include "cusolverSp.h"
#include <cuda_runtime.h>
#include <math.h>

void dpss( size_t N, double NW, double *eigenvector );

int main()
{
    // parameters for generation of dpss
    size_t N = 128;
    double NW = 1;

    double *eigenvector = NULL;
    eigenvector = new double[ N*sizeof( double ) ];
    dpss( N, NW, eigenvector );

    return 0;
}

void dpss( size_t N, double NW, double *eigenvector )
{
    // define matrix T (NxN) 
    double** T = new double*[ N ];
    for(int i = 0; i < N; ++i)
        T[ i ] = new double[ N ];

    // fill in T as function of ( N, W ) 
    // T is a tridiagonal matrix, i. e., it has diagonal, subdiagonal and superdiagonal
    // the others elements are 0
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if( j == i - 1 ) // subdiagonal
                T[ i ][ j ] = ( (double)N - i )*i/2;
            else if( j == i ) // diagonal
                T[ i ][ j ] = pow( (double)(N-1)/2 - i, 2 )*std::cos( 2*NW/(double)N*M_PI )*( j == i );
            else if( j == i + 1 ) // superdiagonal
                T[ i ][ j ] = ( i + 1 )*( (double)N - 1 - i )/2*( j == i + 1 );
            else // others elements
                T[ i ][ j ] = 0;
        }
    }

    // declarations needed
    cusolverStatus_t statCusolver = CUSOLVER_STATUS_SUCCESS;
    cusolverSpHandle_t handleCusolver = NULL;
    cusparseHandle_t handleCusparse = NULL;
    cusparseMatDescr_t descrA = NULL;
    int *h_cooRowIndex = NULL, *h_cooColIndex = NULL;
    double *h_cooVal = NULL; 
    int *d_cooRowIndex = NULL, *d_cooColIndex = NULL, *d_csrRowPtr = NULL; 
    double *d_cooVal = NULL; 
    int nnz; 
    double *h_eigenvector0 = NULL, *d_eigenvector0 = NULL, *d_eigenvector = NULL;
    double max_lambda;

    // define interval of eigenvalues of T
    // interval is [-max_lambda,max_lambda]
    max_lambda = ( N - 1 )*( N + 2 ) + N*( N + 1 )/8 + 0.25;

    // amount of nonzero elements of T
    nnz = 3*N - 2;

    // allocate host memory
    h_cooRowIndex = new int[ nnz*sizeof( int ) ];
    h_cooColIndex = new int[ nnz*sizeof( int ) ];
    h_cooVal = new double[ nnz*sizeof( double ) ];
    h_eigenvector0 = new double[ N*sizeof( double ) ];

    // fill in vectors that describe T as a sparse matrix
    int counter = 0;
    for (int i = 0; i < N; i++ ) {
        for( int j = 0; j < N; j++ ) {
            if( T[ i ][ j ] != 0 ) {
                h_cooColIndex[counter] = j;
                h_cooRowIndex[counter] = i;
                h_cooVal[counter++] = T[ i ][ j ];
            }
        }
    }

    // fill in initial eigenvector guess  
    for( int i = 0; i < N; i++ )
        h_eigenvector0[ i ] =  (double)1/(i+1);

    // allocate device memory
    CUDA_CALL( cudaMalloc((void**)&d_cooRowIndex,nnz*sizeof( int )) ); 
    CUDA_CALL( cudaMalloc((void**)&d_cooColIndex,nnz*sizeof( int )) ); 
    CUDA_CALL( cudaMalloc((void**)&d_cooVal, nnz*sizeof( double )) );
    CUDA_CALL( cudaMalloc((void**)&d_csrRowPtr, (N+1)*sizeof( int )) );
    CUDA_CALL( cudaMalloc((void**)&d_eigenvector0, N*sizeof( double )) );
    CUDA_CALL( cudaMalloc((void**)&d_eigenvector, N*sizeof(d_eigenvector[0])) );

    // copy data to device
    CUDA_CALL( cudaMemcpy( d_cooRowIndex, h_cooRowIndex, (size_t)(nnz*sizeof( int )), cudaMemcpyHostToDevice ) );
    CUDA_CALL( cudaMemcpy( d_cooColIndex, h_cooColIndex, (size_t)(nnz*sizeof( int )), cudaMemcpyHostToDevice ) );
    CUDA_CALL( cudaMemcpy( d_cooVal, h_cooVal, (size_t)(nnz*sizeof( double )), cudaMemcpyHostToDevice ) );
    CUDA_CALL( cudaMemcpy( d_eigenvector0, h_eigenvector0, (size_t)(N*sizeof( double )), cudaMemcpyHostToDevice ) );

    // initialize cusparse and cusolver
    CUSOLVER_CALL( cusolverSpCreate( &handleCusolver ) );
    CUSPARSE_CALL( cusparseCreate( &handleCusparse ) );

    // create and define cusparse matrix descriptor
    CUSPARSE_CALL( cusparseCreateMatDescr(&descrA) );
    CUSPARSE_CALL( cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL ) );
    CUSPARSE_CALL( cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ZERO ) );

    // transform from coordinates (COO) values to compressed row pointers (CSR) values
    CUSPARSE_CALL( cusparseXcoo2csr( handleCusparse, d_cooRowIndex, nnz, N, d_csrRowPtr, CUSPARSE_INDEX_BASE_ZERO ) );

    // define some parameters and call cusolverSpScsreigvsi
    int maxite = 1e6;
    double tol = 1;
    double mu = 0;
    statCusolver = cusolverSpDcsreigvsi( handleCusolver, N ,nnz, descrA, d_cooVal, d_csrRowPtr, d_cooColIndex, max_lambda, d_eigenvector0, maxite, tol, &mu, d_eigenvector );
// here statCusolver = CUSOLVER_INTERNAL_ERROR

    cudaDeviceSynchronize();
    CUDA_CALL( cudaGetLastError() );

    // copy from device to host
    CUDA_CALL( cudaMemcpy( h_eigenvector0, d_eigenvector, (size_t)(N*sizeof( double )), cudaMemcpyDeviceToHost ) );

    // destroy and free stuff
    CUSPARSE_CALL( cusparseDestroyMatDescr( descrA ) );
    CUSPARSE_CALL( cusparseDestroy( handleCusparse ) );
    CUSOLVER_CALL( cusolverSpDestroy( handleCusolver ) );
    CUDA_CALL( cudaFree( d_cooRowIndex ) );
    CUDA_CALL( cudaFree( d_cooColIndex ) );
    CUDA_CALL( cudaFree( d_cooVal ) );
    CUDA_CALL( cudaFree( d_csrRowPtr ) );
    CUDA_CALL( cudaFree( d_eigenvector0 ) );
    CUDA_CALL( cudaFree( d_eigenvector ) );
    delete[] h_eigenvector0;
    delete[] h_cooRowIndex;
    delete[] h_cooColIndex;
    delete[] h_cooVal;
}

Я уже пробовал разные варианты для начального предположения собственного значения (а именно, max_lambda - или mu0 в учебнике по библиотеке cuSolver), начального предположения собственного вектора (h_eigenvector0 или d_eigenvector0), допуска (tol), даже количества максимальной итерации (maxite).

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

Я не знаю, что еще я могу сделать, но если кто-то знает, пожалуйста, дайте мне знать!!

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


person Leonardo Miquelutti    schedule 07.03.2015    source источник


Ответы (1)


Мне кажется (мне), что документация cuSolver может быть неверной в отношении параметра mu.

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

Если я изменю его на указатель устройства:

double mu = 0;
double *d_mu;
CUDA_CALL(cudaMalloc(&d_mu, sizeof(double)));
CUDA_CALL(cudaMemset(d_mu, 0, sizeof(double)));
CUSOLVER_CALL(cusolverSpDcsreigvsi( handleCusolver, N ,nnz, descrA, d_cooVal, d_csrRowPtr, d_cooColIndex, max_lambda, d_eigenvector0, maxite, tol, d_mu, d_eigenvector ));

...
CUDA_CALL(cudaMemcpy(&mu, d_mu, sizeof(double), cudaMemcpyDeviceToHost));

Я могу запустить версию вашего кода без ошибок API или cuda-memcheck ошибок. (За результаты не ручаюсь.)

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

person Robert Crovella    schedule 07.03.2015
comment
Я реализовал около 60 функций cublas и перепроверил все с результатами Matlab, и все мои результаты были достаточно близки к результатам Matlab. Впервые они достаточно разные, особенно собственное значение. Если вы измените начальное предположение собственного вектора на 1/(abs( i - N/2) + 1 ), собственные векторы, которые я получил (изменив N и NW), были близки к векторам Matlab, хотя и не так близко, к которым я привык. Но это сильно зависит от начального предположения собственного вектора (я просто знаю, как должен выглядеть мой собственный вектор в конце, поэтому я принудительно использовал это начальное условие с ранее упомянутым изменением). - person Leonardo Miquelutti; 09.03.2015
comment
Мне жаль. Собственное значение правильное. Собственный вектор, в зависимости от параметров N и NW, может быть самым разным. - person Leonardo Miquelutti; 09.03.2015
comment
Кроме того, я проверил, что результирующий собственный вектор, выдаваемый cusolver, имеет норму, равную 1, как и должно быть. - person Leonardo Miquelutti; 09.03.2015
comment
Если вы хотите предоставить короткий код Matlab, который вы используете для проверки, я посмотрю, хотя я не уверен, что у меня будут какие-либо полезные наблюдения. Я просто пытался продублировать создание вашей T-матрицы в Matlab с такими же вложенными циклами for и с использованием eig - person Robert Crovella; 09.03.2015
comment
Это то, что я только что сделал. Есть еще один способ проверить это с помощью функции dpss, поскольку result = dpss( N, NW, 1 ); - я пытаюсь сгенерировать dpss нулевого порядка, также известного как последовательность Слепяна или нулевого порядка. Однако вполне возможно, что Matlab сгенерирует dpss на основе собственных значений T в соответствии со ссылкой на алгоритм на его веб-странице, что мы и делаем. Любой совет о том, как я могу проверить правильный результат, т.е. д., Matlab или cuSolver? - person Leonardo Miquelutti; 10.03.2015