Hi,

I am currently using MAGMA library to compute the eigen vectors of a symmetric matrix. Unfortunately, the results i got are different than the results of matlab when using "eig" function.

Can you explain why?

This is the code i'm using:

int compute_eigen_vectors (double * covariance_matrix, int dim, double * eigen_matrix,double * covariance_matrix_copy) /

{

cudaMemcpy(covariance_matrix_copy,covariance_matrix,dim*dim*sizeof(double),cudaMemcpyDeviceToDevice);

magma_init (); // initialize Magma

magma_queue_t queue = NULL ;

magma_int_t dev =0;

double * h_work ;

magma_queue_create (dev , & queue );

magma_int_t n=dim , n2=n*n;

magma_int_t lwork ; // h_work size

magma_int_t * iwork ; // workspace

magma_int_t liwork ; // iwork size

double *w1 ; // w1 ,w2 - vectors of eigenvalues

double *r ;

magma_dmalloc_cpu (&w1 ,n);

magma_dmalloc_cpu (&r,n2 );

magma_int_t info ;

magma_int_t incr = 1;

double *d_r ;

magma_dmalloc (& d_r ,n2 );

// Query for workspace sizes

double aux_work [1];

magma_int_t aux_iwork [1];

magma_dmalloc_cpu (&r,n2 ); // host memory for r

magma_dsyevd_gpu ( MagmaVec , MagmaUpper, n , d_r , n , w1 , r , n , aux_work , 1 , aux_iwork , -1 , &info );

lwork = ( magma_int_t ) aux_work [0];

magma_dmalloc_cpu (& h_work , lwork ); // memory for workspace

liwork = aux_iwork [0];

iwork =( magma_int_t *) malloc ( liwork * sizeof ( magma_int_t ));

magma_dmalloc_cpu (& h_work , lwork ); // memory for workspace

lwork = ( magma_int_t ) aux_work [0];

liwork = aux_iwork [0];

iwork = ( magma_int_t *) malloc ( liwork * sizeof ( magma_int_t ));

magma_dmalloc_cpu (& h_work , lwork ); // memory for workspace

magma_dsetmatrix ( n, n, covariance_matrix_copy, n, d_r , n, queue ); // copy a -> d_r

magma_dsyevd_gpu(MagmaVec,MagmaUpper,n,d_r,n,w1,r,n,h_work,lwork,iwork,liwork,&info);

cError = cudaMemcpy(eigen_matrix,d_r,dim*dim*sizeof(double),cudaMemcpyDeviceToDevice);

gpuErrorCheck(cError);

free (w1); // free memory

free (r); // free memory

free ( h_work ); // free memory

magma_free (d_r );

magma_queue_destroy ( queue ); // destroy queue

magma_finalize (); // finalize Magma

return EXIT_SUCCESS ;

}

Thanks

## Coputing eigen vectors

### Re: Coputing eigen vectors

When you say that you get different results than Matlab, can you explain? Can you give a small example, like the input & output for a 4x4 matrix showing the difference? From a cursory view of your code, it looks ok (with caveats below that shouldn't affect the results). So I don't really know what the issue is that you are experiencing.

Note that eigenvectors are not unique. You can multiply by a scalar, or even if normalized, by -1. If there is a multiple eigenvalue, then any linear combination of their eigenvectors is also an eigenvector. The order of eigenvalues may also be different, although it appears that both MAGMA and Matlab will sort them in ascending order. So there are legitimate mathematical reasons why MAGMA and Matlab would produce different — both correct — results. You can check results using formulas that should be around 1e-15:

% assuming symmetric/Hermitian A

>> [m, n] = size( A );

>> [V, D] = eig( A );

>> norm( V*D*V' - A ) / ( norm( A ) * n ) % backward error

ans =

3.5694e-16

>> norm( V' * V - eye(n, n) ) / n % orthogonality error

ans =

2.8968e-16

A few comments on the code:

1) When querying for workspace, I would set both lwork = -1 and liwork = -1, for clarity. You have lwork = 1.

2) You repeat a number of lines which cause memory leaks, for instance:

magma_dmalloc_cpu (&r,n2 ) // occurs 2x

magma_dmalloc_cpu (& h_work , lwork ) // occurs 3x

3) If you use magma_*malloc_cpu, you should pair it with magma_free_cpu instead of free. On some platforms, not pairing them will cause a crash.

Note that eigenvectors are not unique. You can multiply by a scalar, or even if normalized, by -1. If there is a multiple eigenvalue, then any linear combination of their eigenvectors is also an eigenvector. The order of eigenvalues may also be different, although it appears that both MAGMA and Matlab will sort them in ascending order. So there are legitimate mathematical reasons why MAGMA and Matlab would produce different — both correct — results. You can check results using formulas that should be around 1e-15:

% assuming symmetric/Hermitian A

>> [m, n] = size( A );

>> [V, D] = eig( A );

>> norm( V*D*V' - A ) / ( norm( A ) * n ) % backward error

ans =

3.5694e-16

>> norm( V' * V - eye(n, n) ) / n % orthogonality error

ans =

2.8968e-16

A few comments on the code:

1) When querying for workspace, I would set both lwork = -1 and liwork = -1, for clarity. You have lwork = 1.

2) You repeat a number of lines which cause memory leaks, for instance:

magma_dmalloc_cpu (&r,n2 ) // occurs 2x

magma_dmalloc_cpu (& h_work , lwork ) // occurs 3x

3) If you use magma_*malloc_cpu, you should pair it with magma_free_cpu instead of free. On some platforms, not pairing them will cause a crash.

### Re: Coputing eigen vectors

Thanks!

I corrected the code. I can't explain why but it works.

I appreciate your help.

I corrected the code. I can't explain why but it works.

I appreciate your help.