eigenvalue and eigenvector calculation with cusolver from cuda 7.0 RC

Question

I'm trying to calculate the largest eigenvalue/eigenvector pair with cuSolver released in CUDA 7.0 RC. The problem is that I'm getting CUSOLVER_INTERNAL_ERROR, and I don't know what can I do about it.

This is my handy stuff, used to call cuda/cusparse/cusolver functions.

// 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)

And this is my code

#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;
}

I already tried different choices for initial eigenvalue guess (namely max_lambda - or mu0 in the cuSolver library tutorial), initial eigenvector guess (h_eigenvector0 or d_eigenvector0), tolerance (tol), even amount of maximum iteration (maxite).

I already checked if the sparse matrix was properly written (it seemed correct to me). I also checked the returned eigenvector with Matlab, and they are completely different (and I think they shouldn't be).

I don't know what else can I do, but if someone does, please let me konw!!

Thanks in advance.

You're [supposed to provide](http://stackoverflow.com/help/on-topic) an [MCVE](http://stackoverflow.com/help/mcve). — Robert Crovella, Mar 07 '15 at 18:30
I'm sorry guys. My first question asked. I'll do it properly next question I post. Although I found important to post it all due to the function I was having problem is a 13-parameters one, and I imagined that it could be any. — Leonardo Miquelutti, Mar 09 '15 at 00:30

score 3 · Accepted Answer · answered Mar 07 '15 at 22:35

3

It appears (to me) that the cuSolver documentation may be incorrect with respect to the mu parameter.

The documentation appears to indicate that this is in the host memory space, i.e. the 2nd to last parameter should be a host pointer.

If I change it to be a device pointer:

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));

I can get a version of your code to run without any API errors or cuda-memcheck errors. (I can't vouch for the results.)

I have already filed a documentation inquiry with NVIDIA. If you can confirm that you get sensible results with this change, that might also be useful. I took a look at the resultant eigenvalue and eigenvector, and although they did not appear to be bizarre, I could not get them to exactly correlate with my attempt to do the same thing in Octave.

answered Mar 07 '15 at 22:35

Robert Crovella

143,785
11
213
257

I implemented about 60 cublas functions and cross checked all with Matlab results, and all my results were close enough to Matlab ones. For the first time they are reasonably different, especially the eigenvalue value. If you change the initial eigenvector guess to 1/( abs( i - N/2) + 1 ), the eigenvectors I got (changing N and NW) were close to Matlab ones, although not the "close" I'm used to. But it is heavily dependent on the initial eigenvector guess (I just know how my eigenvector should looks like in the end, so I forced this initial condition with the change previously mentioned). – Leonardo Miquelutti Mar 09 '15 at 00:34
I'm sorry. Eigenvalue is correct. Eigenvector, depending on parameters N and NW, can be very different. – Leonardo Miquelutti Mar 09 '15 at 13:25
Also, I checked that the resulting eigenvector outputted by the cusolver has norm equals 1, as it should be. – Leonardo Miquelutti Mar 09 '15 at 19:22
If you want to provide a short matlab code that you use for validation, I'll take a look, although I'm not sure I'll have any useful observations. I was just trying to duplicate your T matrix creation in matlab with the same kind of nested for loops, and using `eig` – Robert Crovella Mar 09 '15 at 20:17
That's what I just did. There is another way of checking it is using dpss function, as `result = dpss( N, NW, 1 );` - I'm trying to generate a dpss of zeroth-order, aka Slepian sequence or zeroth-order. However, it is possible that Matlab generate a dpss based on the eigenvalues of T, accordingly to the algorithm reference on its webpage, which is exactly what we are doing. Any tip on how can I check the correct result, i. e., Matlab or cuSolver? – Leonardo Miquelutti Mar 10 '15 at 20:56

eigenvalue and eigenvector calculation with cusolver from cuda 7.0 RC

1 Answers1