#ifndef maxdim
#define maxdim 130
#endif

#define hugenum 100000000000000000000000.0

extern "C"
{

#ifndef fcomplex_def
#include "fcomplexdef.h"
#endif

void sgeev(float[][maxdim],int*,int*, complex[], 
           complex[][maxdim],int*, float[], int*, int*);
};



class sgeevClass 		// sgeev computes eigensystem for us
{
 float work[2*maxdim];
 double abse[maxdim];
 
 public:
 
 sgeevClass(int dim_);
 int dim;
 void set(int dim_);
 void get_eigensys(float A[][maxdim],int compute_eigvecs=1); 
 void run(float A[][maxdim],int cevc=1) { get_eigensys(A,cevc); }; 
 void get_order();
 int order[maxdim];
 void select(int); 		// selects eigensys
 complex sel_val;
 complex sel_vec[maxdim]; 
 
// sgeev parameters:
 
 int lda;
 int n; 
 complex E[maxdim];			//output
 complex V[maxdim][maxdim];		//output
 int ldv;
 
 int job;     			
 int info;				//output
};

 sgeevClass::sgeevClass(int dim_) { set(dim_); }

void sgeevClass::set(int dim_) 
 { n=dim=dim_; ldv=lda=130; for (int i=0;i<dim;i++) order[i]=i;}
 
 void sgeevClass::get_eigensys(float A[][maxdim],int compute_eigvecs)  
 {
   job=compute_eigvecs;
   sgeev(A,&lda,&n,E,V,&ldv,work,&job,&info);
 }
 
void sgeevClass::get_order() 
{
 int i,j,k,ik,jj;
 double up;
 for (i=0;i<dim;i++) 
 { 
    abse[i]=E[i].im*E[i].im+E[i].re*E[i].re;
    order[i]=-1;
 }
 
 for (i=0;i<dim;i++) 
 {
  for (up=hugenum,k=0;k<dim;k++)
    if (abse[k]<=up) 
     {
        for (jj=1,j=0;j<i;j++) if (order[j]==k) jj=0;  
        if (jj)
        {
           ik=k; 
           up=abse[k];
        }  
     }
  order[i]=ik; 
 }

}
 

void sgeevClass::select(int no) 
{
 int i;
 if (no>=dim) return;
 sel_val.re=E[order[no]].re;
 sel_val.im=E[order[no]].im; 
 for (i=0;i<dim;i++)
 {
   sel_vec[i].re=V[order[no]][i].re; 
   sel_vec[i].im=V[order[no]][i].im; 
 }
}