MAGMA Forum

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

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.

Thanks!

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