/*	Implementation of the LAPACK routine zheev in C++ for
 *	finding eigenvectors and eigenvalues of a square
 *	hermitian matrix H.
 *	
 *	This example contains two versions of the implementation,
 *	not very different from each other, zheeva and zheevb.
 *	
	zheeva assumes the matrix in array form, with elements
	passed C-like, as lines, but complex conjugated. It 
	stores the lower triangle. The trick is that addressing
	matrix elements line-wise is faster in C, so we do it
	like that, just complex conjugating the elements and
	keeping the lower triangle, which amounts to the same
	thing in the end.

	zheevb takes the matrix in 2D form. The elements are
	then  passed Fortran-like as columns to the array.

	in both routines, the initial array is the one used to
	store the eigenvectors later. Thus we spare some time
	with the memory allocation. The eigenvectors are stored,
	in both cases, C-like, as lines.

	sorter is a home-made routine which just finds the largest
	'max' eigenvalues. It's not the brightest, but it does the
	job. It reshuffles the eigenvector matrix accordingly.

	E is the eigenvalue array
	H (or a) is the matrix to be diagonalized, and after
	the diagonalization, the eigenvalue matrix
	max is the number of the largest eigenvalues needed.
	n is the dimension of the matrix to be diagonalized


  Alex Cojuhovschi
  20 March 2006,
  with the due credits to Armand Niederberger

*/

#include <math.h>
#include <complex.h>

void zheeva(complex<double> *H, int n, double *E);
void zheevb(complex<double> **H, int n, double *E);
extern "C" void zheev_(char *jobz, char *uplo, int *n, complex<double> *a,int *lda, double *w, complex<double> *work, int *lwork, double *rwork, int *info);

int mini(int a, int b){
	//computes min(a,b)
	int swapi = a;
	
	if(a > b){
	    swapi = b;
	}
	
	return swapi;
}


void sorter(double *E, complex<double> *H, int max, int range){
	//range is the whole range of eigenvalues
	//max is the max number of eigenvalues we need
	
	int tempint;
	complex<double> ctemp;
	    
	//if the 'max' number is larger than the 'range', then look only for 'range' values
	tempint=mini(range,max);
	
	for(int j = 0; j < tempint; j++){
		int it = j;
		double temp = E[j];
		
		for(int i = j; i < range; i++){
			if(E[i] > temp){
				temp = E[i];
				it = i;
			}
		}
		
		E[it] = E[j];
		E[j] = temp;

		for(int k = 0; k < range; k++){
			ctemp = H[range*j+k];
			H[range*j+k] = H[range*it+k];
			H[range*it+k] = ctemp;
		}
	}
		
		
}    


void zheeva(complex<double> *a, int n, int max, double *E)
{
  char jobz, uplo;
  int info;

	jobz = 'V'; // V/N indicates that eigenvectors should/should not be calculated
	uplo = 'L'; // U/L indicated that the upper/lower triangle of the symmetric matrix is stored
  
	int lwork = 3*n-1;
  
	complex<double> *work = new complex<double>[lwork]; //The work array to be used by zheev an its size
	double *rwork = new double[3*n-2];
  
	zheev_(&jobz, &uplo, &n, a, &n, E, work, &lwork, rwork, &info);//perform the diagonalization

  
	delete [] work;
	delete [] rwork;
  
	sorter(E,a,max,n);

}

void zheevb(complex<double> **H, int n, int max, double *E)
{
  char jobz, uplo;
  int info;

	jobz = 'V'; // V/N indicates that eigenvectors should/should not be calculated
	uplo = 'U'; // U/L indicated that the upper/lower triangle of the symmetric matrix is stored
  
	int lwork = 3*n-1;
  
	complex<double> *work = new complex<double>[lwork]; //The work array to be used by zheev an its size
	double *rwork = new double[3*n-2];
	complex<double> *a = new complex<double>[n*n];
	
	for (int j=0; j < n; j++){
		for (int i=0; i < n; i++){
			a[n*j+i] = H[i][j];//passed as rods
		}
	}
	
  
	zheev_(&jobz, &uplo, &n, a, &n, E, work, &lwork, rwork, &info);//perform the diagonalization

  
	delete [] work;
	delete [] rwork;

	sorter(E,a,max,n);//Sort them
	
	for (int i=0; i < n; i++){//Pass them back to the eigenvectors matrix
		for(int j = 0; j < n; j++){
			H[i][j] = a[n*i+j];//passed as lines
		}
	}
  
	delete [] a;


}