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