Thank you very much Dr Klaus for your reply.
Attached my c code which i try to calculate the eigenvector of two matrix
AA and BB. For the eigenvalue it is ok but for the eigenvector it is
completely different from matlab.
Thank you for your help.
Best regards
2014-05-31 1:26 GMT+03:00 Klaus Huthmacher <[email protected]>:
> Dear Belwafi Kais,
>
> maybe you can give us some sourcecode which we can check, e.g. the matrix
> you are working with?
>
> Furthermore, from my humble experience, some eigensolve programs give back
> normalized eigenvectors some not.
>
> Regards,
> -- Klaus.
>
>
--
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
void plot_matrix_complex__(gsl_matrix_complex *Matrix__,int t1)
{
int i,j;
for (i = 0; i < t1; i++)
{
gsl_vector_complex_view evec_i= gsl_matrix_complex_column (Matrix__, i);
for (j = 0; j < t1; ++j)
{
gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
printf("%.2g \t", GSL_REAL(z));
}
printf (";\n");
}
}
void plot_matrix__(gsl_matrix *Data_matrix, int t1)
{
int i,j;
for (i = 0; i < t1; i++)
{
for (j = 0; j < t1; ++j)
printf("%.2g \t", gsl_matrix_get(Data_matrix,i,j));
printf (";\n");
}
}
int main (void)
{
int n=4;
gsl_matrix *trial1 = gsl_matrix_alloc (n,n);
gsl_matrix *trial2 = gsl_matrix_alloc (n,n);
gsl_eigen_genv_workspace * Trials_1_2 = gsl_eigen_genv_alloc (n);
gsl_vector_complex *alpha = gsl_vector_complex_alloc (n);
gsl_vector *beta= gsl_vector_alloc (n);
gsl_matrix_complex *eigvector=gsl_matrix_complex_alloc (n,n);
gsl_matrix *Q = gsl_matrix_alloc (n,n);
gsl_matrix *Z = gsl_matrix_alloc (n,n);
// Initiliaze data
{
int i,j;
for (i = 0; i < trial1->size1; i++)
for (j = 0; j < trial1->size2; j++)
{
gsl_matrix_set (trial1, i, j,1.5 );//0.23 + 1*i + j
gsl_matrix_set (trial2, i, j,3.26);//0.50 + 7*j
}
double AA[16] = { 8, -5, -5, -7, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10 };
double BB[16] = { 7, -10, -9, -1, 4, 3, -7, 9, -4, 1, -9, 9, -7, 5, -5, -1} ;
gsl_matrix_view A = gsl_matrix_view_array(AA, 4, 4);
gsl_matrix_view B = gsl_matrix_view_array(BB, 4, 4);
gsl_eigen_genv(&A.matrix,&B.matrix,alpha,beta,eigvector,Trials_1_2);
for (j = 0; j < trial1->size1; j++)
{
//gsl_complex eval_i= gsl_vector_complex_get (alpha, j);
if(gsl_vector_get (beta, j)!=0)
printf ("eigenval(%d) = %g\n",j,GSL_REAL(gsl_vector_complex_get (alpha, j))/gsl_vector_get (beta, j));
else
if(GSL_REAL(gsl_vector_complex_get (alpha, j))!=0)
printf ("eigenval(%d) = %g\n",j,gsl_vector_get (beta, j))/GSL_REAL(gsl_vector_complex_get (alpha, j));
else
printf ("eigenval(%d) = %g\n",j,GSL_REAL(gsl_vector_complex_get (alpha, j))/gsl_vector_get (beta, j));
// printf ("eigenval(%d) = %g\n",j,GSL_REAL(gsl_vector_complex_get (alpha, j)));
}
printf ("eigenvector = \n");
plot_matrix_complex__(eigvector,trial1->size1);
printf ("-----------------------------------------------\n");
printf ("**************** Q ***************\n");
gsl_eigen_genv_QZ(trial1,trial2,alpha,beta,eigvector,Q,Z,Trials_1_2);
plot_matrix__(Q,trial1->size1);
printf ("**************** Z ***************\n");
plot_matrix__(Z,trial1->size1);
printf ("______________ eigenvector __________\n");
plot_matrix_complex__(eigvector,trial1->size1);
printf ("-----------------------------------------------\n");
}
gsl_matrix_free(Q);
gsl_matrix_free(Z);
gsl_matrix_free(trial1);
gsl_matrix_free(trial2);
gsl_vector_complex_free(alpha);
gsl_vector_free(beta);
gsl_eigen_genv_free(Trials_1_2 );
return 0;
}
