/***********************************************/
/* T.Kouya's GSL sample program collection     */
/*      Eingenvalue Problem                    */
/*                   for Real Symmetrix Matrix */
/*                                             */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2007-10-04                    */
/***********************************************/
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sort_vector.h>

/* Dimension of Matrix and Vectors */
#define DIM 5

int main(void)
{
	int i, j;
	gsl_permutation *perm;
	gsl_matrix *frank_a, *matrix, *eig_vecs, *tmp_matrix;
	gsl_vector *eig_vec, *eig;
	gsl_eigen_symm_workspace *workspace;
	gsl_eigen_symmv_workspace *workspace_v;

	/* allocate a, x, b */
	frank_a = gsl_matrix_alloc(DIM, DIM);
	matrix = gsl_matrix_alloc(DIM, DIM);
	eig_vec = gsl_vector_alloc(DIM);
	eig = gsl_vector_alloc(DIM);

	/* set Frank matrix       */
	/* a_ij = DIM - min(i,j) + 1 */
	for(i = 0; i < DIM; i++)
		for(j = 0; j < DIM; j++)
			gsl_matrix_set(frank_a, i, j, (double)DIM - (double)GSL_MAX(i, j) );

	/* Print matrix */
	printf("Frank Matrix(DIM = %d)\n", DIM);
	for(i = 0; i < DIM; i++)
	{
		printf("%3d: ", i);
		for(j = 0; j < DIM; j++)
			printf("%g ", gsl_matrix_get(frank_a, i, j));
		printf("\n");
	}
	printf("\n");

/* Get Eigensystem of Frank Matrix : Eigenvalues ONLY */

	/* Initialize workspace and copy Frank matrix to "matrix" */
	workspace = gsl_eigen_symm_alloc(DIM);
	gsl_matrix_memcpy(matrix, frank_a);

	/* get eigenvalues (matrix is destroyed) */
	gsl_eigen_symm(matrix, eig, workspace);

	/* sort eigenvalues */
	gsl_sort_vector(eig);

	/* print eigenvalues */
	printf(" i       eigenvalues        \n");
	for(i = 0; i < DIM; i++)
		printf("%3d %25.17e\n", i, gsl_vector_get(eig, i));

	/* free workspace */
	gsl_eigen_symm_free(workspace);

/* Get Eigensystem of Frank Matrix : Eigenvalues and Eigenvectors */

	/* Initialize workspace and eigenvectors stored, and set matrix */
	eig_vecs = gsl_matrix_alloc(DIM, DIM);
	workspace_v = gsl_eigen_symmv_alloc(DIM);
	gsl_matrix_memcpy(matrix, frank_a);

	/* get eigenvalues and its eigenvectors (matrix is destroyed) */
	gsl_eigen_symmv(matrix, eig, eig_vecs, workspace_v);

	/* sort eigenvalues with permutation */
	perm = gsl_permutation_alloc(DIM);
	gsl_sort_vector_index(perm, eig);

	/* print eigenvalues */
	printf(" i       eigenvalues        \n");
	for(i = 0; i < DIM; i++)
		printf("%3d %25.17e\n", i, gsl_vector_get(eig, gsl_permutation_get(perm, i)));
	printf("\n");

	/* print eigenvectors */
	printf("eigenvectors : [v1 v2 ... vn]\n");
	for(i = 0; i < DIM; i++)
		printf("eig%2d = %5.1e  ", i, gsl_vector_get(eig, i));
	printf("\n");
	for(i = 0; i < DIM; i++)
	{
		for(j = 0; j < DIM; j++)
			printf("%16.9e ", gsl_matrix_get(eig_vecs, i, j));
		printf("\n");
	}
	printf("\n");

	/* free permutation */
	gsl_permutation_free(perm);

	/* A * [v1 v2 ... vn] - diag[eig1 eig2 ... eign] * [v1 v2 ... vn] == 0 ? */
	tmp_matrix = gsl_matrix_alloc(DIM, DIM);
	gsl_matrix_set_identity(tmp_matrix);
	for(i = 0; i < DIM; i++)
		for(j = 0; j < DIM; j++)
			gsl_matrix_set(matrix, j, i, gsl_vector_get(eig, i));

	gsl_matrix_mul_elements(matrix, eig_vecs);
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, frank_a, eig_vecs, -1.0, matrix);

	/* print residual */
	printf("A * [v1 v2 ... vn] - diag[eig1 eig2 ... eign] * [v1 v2 ... vn]\n");
	for(i = 0; i < DIM; i++)
	{
		printf("%3d: ", i);
		for(j = 0; j < DIM; j++)
			printf("%g ", gsl_matrix_get(matrix, i, j));
		printf("\n");
	}
	printf("\n");

	/* free workspace */
	gsl_matrix_free(eig_vecs);
	gsl_matrix_free(tmp_matrix);
	gsl_eigen_symmv_free(workspace_v);

	/* free a, x, b */
	gsl_matrix_free(frank_a);
	gsl_matrix_free(matrix);
	gsl_vector_free(eig_vec);
	gsl_vector_free(eig);

	return 0;
}