/* calculate eigenvalues and eigenvectors of a Hermitian matrix */

#include <math.h>
#include "cmplx.h"
#include "../iolib/message.h"

#define SMALL 1.e-04

float **farray2();
float *farray1();

void eigherm(H,n1,n2,lambda)
complex **H;
int n1,n2;
float *lambda;
{
float **a,*d,*e;
int i,j,n;
double err;
  n=n2-n1+1;
  a=farray2(1,2*n,1,2*n); /* RECIPES package */
  d=farray1(1,2*n);
  e=farray1(1,2*n);
  
  for(i=0;i<n;i++)
  for(j=0;j<n;j++) 
  { a[i+1+n][j+1+n]=a[i+1][j+1]=H[i+n1][j+n1].real;
    a[i+1][j+1+n]=-H[i+n1][j+n1].imag;
    a[i+1+n][j+1]= H[i+n1][j+n1].imag;
  }
  
  tred2(a,2*n,d,e); /* RECIPES package */
  tqli(d,e,2*n,a);  /* RECIPES package */
  
  for(i=0;i<n;i++)
  { lambda[i]=d[2*i+1];
    err=(d[2*i+1]-d[2*i+2])/d[2*i+2];
    if(fabs(err)>SMALL) STOP(eigherm:error 1);
    for(j=0;j<n;j++)
    { H[j+n1][i+n1].real=a[j+1][2*i+1];
      H[j+n1][i+n1].imag=a[j+1+n][2*i+1];
    }
  }
  free_farray2(a,1,2*n,1);
  free_farray1(d,1);
  free_farray1(e,1);
}