/*
 * eigen.c - calculation of eigen values and vectors.
 *
 * (C) Copyright 2001 by NetGroup A/S. All rights reserved.
 *
 * $Log: eigen.c,v $
 * Revision 1.1  2006/06/20 15:57:22  djburke
 * Hopefully a saner way to build Basic/MatrixOps
 *
 * Revision 1.1  2005/01/08 09:22:57  zowie
 * Added non-symmetric matrices to eigens; updated version to 2.4.2cvs.
 *
 * Revision 1.1.1.1  2001/07/06 13:39:35  kneth
 * Initial import of code.
 *
 *
 * Eigen is a library for computing eigenvalues and eigenvectors of general
 * matrices. There is only one routine exported, namely Eigen.
 *
 * The meaning of the arguments to Eigen is:
 *   1.   The dimension of the general matrix (n).
 *   2.   A general matrix (A).   
 *   3.   The maximal number of iterations.
 *   4.   The precision.
 *   5.   A vector with the eigenvalues.
 *   6.   A matrix with the eigenvectors.
 *
 */

#include "complex.h"
#include "matrix.h"
#include <math.h>
#include <stdio.h>

void BlockCheck(double **A, int n, int i, int *block, double epsx) {

  /* block == 1 <=> TRUE, block == 0 <=> FALSE */

  if (i==n)
    *block=0;
  else {
    if ((fabs(A[i-1][i]-A[i][i-1])>epsx) && 
	(fabs(A[i-1][i-1]-A[i][i])<=epsx))
      *block=1;
    else
      *block=0;
  } /* else */
} /* BlockCheck */


void PrintEigen(int n, double **A, double **B, double eps, FILE *outfile) {

  int     i, j;
  int     block;

  fprintf(outfile, "\nEigenvalues:\t\t\tRe\t\t\tIm\n");
  i=1;
  do {
    BlockCheck(A, n, i, &block, eps);
    if (block==1) {
      fprintf(outfile, "\t\t\t\t%e\t\t%e\n", A[i-1][i-1], A[i-1][i]);
      fprintf(outfile, "\t\t\t\t%e\t\t%e\n", A[i][i], A[i][i-1]);
      i+=2;
    } else {
      fprintf(outfile, "\t\t\t\t%e\t\t%e\n", A[i-1][i-1], 0.0);
      i++;
    } /* if else */
  } while (i!=(n+1));
  fprintf(outfile, "\nEigenvectors:\t\t\tRe\t\t\tIm\n");
  i=1;
  do {
    BlockCheck(A, n, i, &block, eps);
    if (block==1) {
      for(j=1; j<=n; j++) 
	fprintf(outfile, "\t\t\t\t%e\t\t%e\n", B[j-1][i-1], B[j-1][i]);
      fprintf(outfile, "\n");
      for(j=1; j<=n; j++)
	fprintf(outfile, "\t\t\t\t%e\t\t%e\n", B[j-1][i-1], -B[j-1][i]);
      fprintf(outfile, "\n");
      i+=2;
    } else {
      for(j=1; j<=n; j++) 
	fprintf(outfile, "\t\t\t\t%e\t\t%e\n", B[j-1][i-1], 0.0);
      fprintf(outfile, "\n");
      i++;
    } /* if else */
  } while (i!=(n+1));
} /* PrintEigen */


void NormalizingMatrix(int n, double **A, int fixedref, int *ref, 
		       double **V, double eps) {

  int      j, col, block;
  SSL_Complex  c1, c2, c3;
  double   cd1, cd2, sqrnorm, norm, normi, max;

  col=1;
  do {
    if (fixedref==0) {
      *ref=1;
      SSL_ComplexAssign(V[*ref-1][col-1], V[*ref-1][col], &c1);
      max=SSL_ComplexNorm(c1);
      for(j=2; j<=n; j++) {
	SSL_ComplexAssign(V[j-1][col-1], V[j-1][col], &c2);
	sqrnorm=SSL_ComplexNorm(c2);
	if (sqrnorm>max) {
	  *ref=j;
	  max=sqrnorm;
	} /* if */
      } /* for j */
    } /* if fixedref */
    BlockCheck(A, n, col, &block, eps);
    if (block==1) {
      SSL_ComplexAssign(V[*ref-1][col-1], V[*ref-1][col], &c1);
      for(j=1; j<=n; j++) {
	SSL_ComplexAssign(V[j-1][col-1], V[j-1][col], &c2);
	SSL_ComplexDiv(c2, c1, &c3);
	V[j-1][col-1]=c3.re;
	V[j-1][col]=c3.im;
      } /* for j */
      col+=2;
    } /* if */
    else {
      norm=fabs(V[*ref-1][col-1]);
      if (norm!=0.0)
	for(j=1; j<=n; j++)
	  V[j-1][col-1]/=norm;
      col++;
    } /* else */
  } while (col<=n);
} /* NormalizingMatrix */

void Permutation(int n, double **P, double **A, double **B, int colon,
		 double eps) {

  int      *nr;
  int      block, OK;
  double   max, y, x;
  int      im, j, ki, u, v, i, k, ii;
  double   **AA;

  nr=IntVectorAlloc(n);
  AA=MatrixAlloc(n);

  MatrixCopy(n, AA, A);
  for(i=1; i<=n; i++) {
    nr[i-1]=i;
    for(k=1; k<=n; k++)
      P[i-1][k-1]=0.0;
  } /* for i */
  i=ii=ki=1;
  while (i<n) {
    BlockCheck(A, n, i, &block, eps);
    if (block==1) {
      A[i][i]=A[i-1][i-1];
      AA[i][i]=AA[i-1][i-1];
      if (A[i-1][i]>0.0) {
	A[i][i-1]=A[i-1][i];
	A[i-1][i]=-A[i][i-1];
	AA[i][i-1]=AA[i-1][i];
	AA[i-1][i]=-AA[i][i-1];
	for(j=1; j<=n; j++) 
	  B[j-1][i]=-B[j-1][i];
      } else {
	A[i][i-1]=-A[i-1][i];
	AA[i][i-1]=-AA[i-1][i];
      } /* else */
      j=i;
      for(k=ii; k<=(ii+1); k++) {
	x=AA[k-1][k-1];
	AA[k-1][k-1]=A[j-1][j-1];
	AA[j-1][j-1]=x;
	u=nr[k-1];
	nr[k-1]=nr[j-1];
	nr[j-1]=u;
	j++;
      } /* for k */
      if (ii>1) {
	if (AA[ii-1][ii-1]>AA[0][0]) {
	  j=ii;
	  for(k=1; k<=2; k++) {
	    x=AA[k-1][k-1];
	    AA[k-1][k-1]=A[j-1][j-1];
	    AA[j-1][j-1]=x;
	    u=nr[k-1];
	    nr[k-1]=nr[j-1];
	    nr[j-1]=u;
	    j++;
	  } /* for k */
	} /* if */
      } /* if */
      ki=i;
      i+=2;
      ii+=2;
    } /* if */
    else
      i++;
  } /* while */

  if (n>3) {
    do {
      im=ii;
      i=ii;
      max=AA[im-1][im-1];
      do {
	i++;
	if (AA[i-1][i-1]>max) {
	  im=i;
	  max=AA[i-1][i-1];
	} /* if */
      } while (i<n);
      if (im>ii) {
	x=AA[ii-1][ii-1];
	u=nr[ii-1];
	AA[ii-1][ii-1]=max;
	nr[ii-1]=nr[im-1];
	AA[im-1][im-1]=x;
	nr[im-1]=u;
      } /* if */
      ii++;
    } while (ii<n);
  } /* if */
  for(i=1; i<=n; i++) {
    if (colon==1)
      P[nr[i-1]-1][i-1]=1.0;
    else
      P[i-1][nr[i-1]-1]=1.0;
  } /* for i */

  MatrixFree(n, AA);
  IntVectorFree(n, nr);
} /* Permutation */


void Swap(int n, double **A, double **B, double epsx) {

  double **PR, **PS;
  double **temp;
  int      i, j;
  
  PR=MatrixAlloc(n);
  PS=MatrixAlloc(n);
  temp=MatrixAlloc(n);

  Permutation(n, PS, A, B, 1, epsx);
  MatrixMul(n, temp, B, PS);
  MatrixCopy(n, B, temp);
  Transpose(n, PR, PS);
  MatrixMul(n, temp, PR, A);
  MatrixCopy(n, A, temp);
  MatrixMul(n, temp, A, PS);
  MatrixCopy(n, A, temp);

  MatrixFree(n, PR);
  MatrixFree(n, PS);
  MatrixFree(n, temp);
} /* Swap */

void Balance(int n, int b, double **a, int *low, int *hi, double *d) {

  int     i, j, k, l;
  double  b2, c, f, g, r, s;
  int     noconv;

  b2=b*b;
  l=1;
  k=n;

  L110:
  for(j=k; j>=1; j--) {
    r=0.0;
    for(i=1; i<=(j-1); i++) 
      r+=fabs(a[j-1][i-1]);
    for(i=(j+1); i<=k; i++)
      r+=fabs(a[j-1][i-1]);
    if (r==0.0) {
      d[k-1]=(double)j; 
      if (j!=k) { 
	for(i=1; i<=k; i++) { 
	  f=a[i-1][j-1]; 
	  a[i-1][j-1]=a[i-1][k-1]; 
	  a[i-1][k-1]=f; 
	} 
	for(i=l; i<=n; i++) { 
	  f=a[j-1][i-1]; 
	  a[j-1][i-1]=a[k-1][i-1]; 
	  a[k-1][i-1]=f; 
	}
      }
      k--;
      goto L110;
    } /* if */
  } /* for j */

  L120:
  for(j=l; j<=k; j++) {
    c=0.0;
    for (i=l; i<=(j-1); i++)
      c+=fabs(a[i-1][j-1]);
    for(i=(j+1); i<=k; i++)
      c+=fabs(a[i-1][j-1]);
    if (c==0.0) {
      d[l-1]=(double)j; 
      if (j!=l) { 
	for(i=1; i<=k; i++) { 
	  f=a[i-1][j-1]; 
	  a[i-1][j-1]=a[i-1][l-1]; 
	  a[i-1][l-1]=f; 
	} 
	for(i=l; i<=n; i++) { 
	  f=a[j-1][i-1]; 
	  a[j-1][i-1]=a[l-1][i-1]; 
	  a[l-1][i-1]=f; 
	} 
      }
      l++;
      goto L120;
    } /* if */
  } /* for j */

  *low=l;
  *hi=k;
  for(i=l; i<=k; i++)
    d[i-1]=1.0;

  L130:
  noconv=0;
  for(i=l; i<=k; i++) {
    r=c=0.0;
    for(j=l; j<=(i-1); j++) {
      c+=fabs(a[j-1][i-1]);
      r+=fabs(a[i-1][j-1]);
    } /* for j */
    for(j=(i+1); j<=k; j++) {
      c+=fabs(a[j-1][i-1]);
      r+=fabs(a[i-1][j-1]);
    } /* for j */
    g=r/((double) b);
    f=1.0;
    s=c+r;

    L140:
    if (c<g) {
      f*=(double) b;
      c*=(double) b2;
      goto L140;
    } /* if */
    g=r*((double) b);
    
    L150:
    if (c>=g) {
      f/=(double) b;
      c/=(double) b2;
      goto L150;
    } /* if */

    if ((c+r)/f<(0.95*s)) {
      g=1.0/f;
      d[i-1]*=f;
      noconv=1;
      for(j=l; j<=n; j++) 
	a[i-1][j-1]*=g;
      for(j=1; j<=k; j++)
	a[j-1][i-1]*=f;
    } /* if */
  } /* for i */
  if (noconv==1)
    goto L130;
} /* Balance */

void BalBak(int n, int low, int hi, int m, double **z, double *d) {

  int     i, j, k;
  double  s;

  for(i=low; i<=hi; i++) {
    s=d[i-1];
    for(j=1; j<=m; j++) 
      z[i-1][j-1]*=s;
  } /* for i */
  for(i=(low-1); i>=1; i--) {
    k=(int)floor(d[i-1]+0.5);
    if (k!=i)
      for(j=1; j<=m; j++) {
	s=z[i-1][j-1];
	z[i-1][j-1]=z[k-1][j-1];
	z[k-1][j-1]=s;
      } /* for j */
  } /* for i */
  for(i=(hi+1); i<=n; i++) {
    k=(int)floor(d[i-1]+0.5);
    if (k!=i)
      for(j=1; j<=m; j++) {
	s=z[i-1][j-1];
	z[i-1][j-1]=z[k-1][j-1];
	z[k-1][j-1]=s;
      } /* for j */
  } /* for i */
} /* BalBak */

void Elmhes(int n, int k, int l, double **a, int *index) {

  int    i, j, la, m;
  double x, y;

  la=l-1;
  for(m=(k+1); m<=la; m++) {
    i=m;
    x=0.0;
    for(j=m; j<=l; j++) 
      if (fabs(a[j-1][m-2])>fabs(x)) {
	x=a[j-1][m-2];
	i=j;
      } /* if */
      index[m-1]=i;
      if (i!=m) {
	for(j=(m-1); j<=n; j++) {
	  y=a[i-1][j-1];
	  a[i-1][j-1]=a[m-1][j-1];
	  a[m-1][j-1]=y;
	} /* for j */
	for(j=1; j<=l; j++) {
	  y=a[j-1][i-1];
	  a[j-1][i-1]=a[j-1][m-1];
	  a[j-1][m-1]=y;
	} /* for j */
      } /* if */
      if (x!=0.0) 
	for(i=(m+1); i<=l; i++) {
	  y=a[i-1][m-2];
	  if (y!=0.0) {
	    a[i-1][m-2]=y/x;
	    y/=x;
	    for(j=m; j<=n; j++)
	      a[i-1][j-1]-=y*a[m-1][j-1];
	    for(j=1; j<=l; j++)
	      a[j-1][m-1]+=y*a[j-1][i-1];
	  } /* if */
	} /* for i */
  } /* for m */
} /* Elmhes */

void Elmtrans(int n, int low, int upp, double **h, int *index, 
	      double **v) {

  int   i, j, k;

  for(i=1; i<=n; i++) {
    for(j=1; j<=n; j++)
      v[i-1][j-1]=0.0;
    v[i-1][i-1]=1.0;
  } /* for i */
  for(i=(upp-1); i>=(low+1); i--) {
    j=index[i-1];
    for(k=(i+1); k<=upp; k++)
      v[k-1][i-1]=h[k-1][i-2];
    if (i!=j) {
      for(k=i; k<=upp; k++) {
	v[i-1][k-1]=v[j-1][k-1];
	v[j-1][k-1]=0.0;
      } /* for k */
      v[j-1][i-1]=1.0;
    } /* if */
  } /* for i */
} /* Elmtrans */


void hqr2(int n, int low, int upp, int maxits, double macheps,
          double **h, double **vecs, double *wr,
          double *wi, int *cnt, int *fail) {

  int     i, j, k, l, m, na, its, en, dummy;
  double  p, q, r, s, t, w, x, y, z, ra, sa, vr, vi, norm;
  int     notlast;
  SSL_Complex c1, c2, c3;

  *fail=0;
  for(i=1; i<=(low-1); i++) {
    wr[i-1]=h[i-1][i-1];
    wi[i-1]=0.0;
    cnt[i-1]=0;
  } /* for i */
  for(i=(upp+1); i<=n; i++) {
    wr[i-1]=h[i-1][i-1];
    wi[i-1]=0.0;
    cnt[i-1]=0;
  } /* for i */
  en=upp;
  t=0.0;

  L210:
  if (en<low)
    goto L260;
  its=0;
  na=en-1;

  L220:
  for(l=en; l>=(low+1); l--)
    if (fabs(h[l-1][l-2])<=
	macheps*(fabs(h[l-2][l-2])+fabs(h[l-1][l-1])))
      goto L231;
  l=low;

  L231:
  x=h[en-1][en-1];
  if (l==en)
    goto L240;
  y=h[na-1][na-1];
  w=h[en-1][na-1]*h[na-1][en-1];
  if (l==na)
    goto L250;
  if (its==maxits) {
    cnt[en-1]=maxits+1;
    *fail=1;
    goto L270;
  } /* if */
  if ((its % 10)==0) {
    t+=x;
    for(i=low; i<=en; i++)
      h[i-1][i-1]-=x;
    s=fabs(h[en-1][na-1])+fabs(h[na-1][en-3]);
    y=0.75*s;
    x=y;
    w=-0.4375*s*s;
  } /* if */
  its++;

  for(m=(en-2); m>=l; m--) {
    z=h[m-1][m-1];
    r=x-z;
    s=y-z;
    p=(r*s-w)/h[m][m-1]+h[m-1][m];
    q=h[m][m]-z-r-s;
    r=h[m+1][m];
    s=fabs(p)+fabs(q)+fabs(r);
    p/=s;
    q/=s;
    r/=s;
    if (m==1)
      goto L232;
    if ((fabs(h[m-1][m-2])*(fabs(q)+fabs(r)))<=
	(macheps*fabs(p)*(fabs(h[m-2][m-2])+fabs(z)+fabs(h[m][m]))))
      goto L232;
  } /* for m */

  L232:
  for(i=(m+2); i<=en; i++)
    h[i-1][i-3]=0.0;
  for(i=(m+3); i<=en; i++)
    h[i-1][i-4]=0.0;

  for(k=m; k<=na; k++) {
    if (k!=na)
      notlast=1;
    else
      notlast=0;
    if (k!=m) {
      p=h[k-1][k-2];
      q=h[k][k-2];
      if (notlast==1) 
	r=h[k+1][k-2];
      else
	r=0.0;
      x=fabs(p)+fabs(q)+fabs(r);
      if (x==0.0) 
	goto L233;
      p/=x;
      q/=x;
      r/=x;
    } /* if */
    s=sqrt(p*p+q*q+r*r);
    if (p<0)
      s=-s;
    if (k!=m)
      h[k-1][k-2]=-s*x;
    else
      if (l!=m)
	h[k-1][k-2]=-h[k-1][k-2];
    p+=s;
    x=p/s;
    y=q/s;
    z=r/s;
    q/=p;
    r/=p;

    for(j=k; j<=n; j++) {
      p=h[k-1][j-1]+q*h[k][j-1];
      if (notlast==1) {
	p+=r*h[k+1][j-1];
	h[k][j-1]-=p*z;
      } /* if */
      h[k][j-1]-=p*y;
      h[k-1][j-1]-=p*x;
    } /* for j */
    if ((k+3)<en)
      j=k+3;
    else
      j=en;
    
    for(i=1; i<=j; i++) {
      p=x*h[i-1][k-1]+y*h[i-1][k];
      if (notlast==1) {
	p+=z*h[i-1][k+1];
	h[i-1][k+1]-=p*r;
      } /* if */
      h[i-1][k]-=p*q;
      h[i-1][k-1]-=p;
    } /* for i */

    for(i=low; i<=upp; i++) {
      p=x*vecs[i-1][k-1]+y*vecs[i-1][k];
      if (notlast==1) {
	p+=z*vecs[i-1][k+1];
	vecs[i-1][k+1]-=p*r;
      } /* if */
      vecs[i-1][k]-=p*q;
      vecs[i-1][k-1]-=p;
    } /* for i */

  L233:
    dummy=0; 
  } /* for k */
  goto L220;
  
  L240:
  h[en-1][en-1]=x+t;
  wr[en-1]=h[en-1][en-1];
  wi[en-1]=0.0;
  cnt[en-1]=its;
  en=na;
  goto L210;

  L250:
  p=0.5*(y-x);
  q=p*p+w;
  z=sqrt(fabs(q));
  h[en-1][en-1]=x+t;
  x=h[en-1][en-1];
  h[na-1][na-1]=y+t;
  cnt[en-1]=-its;
  cnt[na-1]=its;
  if (q>0.0) {
    if (p<0.0)
      z=p-z;
    else
      z+=p;
    wr[na-1]=x+z;
    s=x-w/z;
    wr[en-1]=s;
    wi[na-1]=0.0;
    wi[en-1]=0.0;
    x=h[en-1][na-1];
    r=sqrt(x*x+z*z);
    p=x/r;
    q=z/r;
    for(j=na; j<=n; j++) {
      z=h[na-1][j-1];
      h[na-1][j-1]=q*z+p*h[en-1][j-1];
      
      /* h[en-1][j-1]=q*h[en-1][j-1]-p*z */
      h[en-1][j-1]*=q; 
      h[en-1][j-1]-=p*z; 
    } /* for j */
    for(i=1; i<=en; i++) {
      z=h[i-1][na-1];
      h[i-1][na-1]=q*z+p*h[i-1][en-1];

      /* h[i-1][en-1]=q*h[i-1][en-1]-p*z */
      h[i-1][en-1]*=q; 
      h[i-1][en-1]-=p*z; 
    } /* for i */
    for(i=low; i<=upp; i++) {
      z=vecs[i-1][na-1];
      vecs[i-1][na-1]=q*z+p*vecs[i-1][en-1];

      /* vecs[i-1][en-1]=q*vecs[i-1][en-1]-p*z */
      vecs[i-1][en-1]*=q; 
      vecs[i-1][en-1]-=p*z; 
    } /* for i */
  } /* if */
  else {
    wr[na-1]=x+p;
    wr[en-1]=x+p;
    wi[na-1]=z;
    wi[en-1]=-z;
  } /* else */
  en-=2;
  goto L210;

  L260:
  norm=0.0;
  k=1;
  for(i=1; i<=n; i++) {
    for(j=k; j<=n; j++)
      norm+=fabs(h[i-1][j-1]);
    k=i;
  } /* for i */

  for(en=n; en>=1; en--) {
    p=wr[en-1];
    q=wi[en-1];
    na=en-1;
    if (q==0.0) {
      m=en;
      h[en-1][en-1]=1.0;
      for(i=na; i>=1; i--) {
	w=h[i-1][i-1]-p;
	r=h[i-1][en-1];
	for(j=m; j<=na; j++)
	  r+=h[i-1][j-1]*h[j-1][en-1];
	if (wi[i-1]<0.0) {
	  z=w;
	  s=r;
	} /* if */
	else {
	  m=i;
	  if (wi[i-1]==0.0)  {
	    if (w!=0.0)
	      h[i-1][en-1]=-r/w;
	    else
	      h[i-1][en-1]=-r/macheps/norm;
	  } else {
	    x=h[i-1][i];
	    y=h[i][i-1];
	    q=pow(wr[i-1]-p, 2.0)+wi[i-1]*wi[i-1];
	    t=(x*s-z*r)/q;
	    h[i-1][en-1]=t;
	    if (fabs(x)>fabs(z))
	      h[i][en-1]=(-r-w*t)/x;
	    else
	      h[i][en-1]=(-s-y*t)/z;
	  } /* else */
	} /* else */
      } /* i */
    } else
      if (q<0.0) {
	m=na;
	if (fabs(h[en-1][na-1])>fabs(h[na-1][en-1])) {
	  h[na-1][na-1]=-(h[en-1][en-1]-p)/h[en-1][na-1];
	  h[na-1][en-1]=-q/h[en-1][na-1];
	} /* if */ 
	else {
	  SSL_ComplexAssign(-h[na-1][en-1], 0.0, &c1);
	  SSL_ComplexAssign(h[na-1][na-1]-p, q, &c2);
	  SSL_ComplexDiv(c1, c2, &c3);
	  h[na-1][na-1]=c3.re;
	  h[na-1][en-1]=c3.im;
	} /* else */
	h[en-1][na-1]=1.0;
	h[en-1][en-1]=0.0;
	for(i=(na-1); i>=1; i--) {
	  w=h[i-1][i-1]-p;
	  ra=h[i-1][en-1];
	  sa=0.0;
	  for(j=m; j<=na; j++) {
	    ra+=h[i-1][j-1]*h[j-1][na-1];
	    sa+=h[i-1][j-1]*h[j-1][en-1];
	  } /* for j */
	  if (wi[i-1]<0.0) {
	    z=w;
	    r=ra;
	    s=sa;
	  } /* if */
	  else {
	    m=i;
	    if (wi[i-1]==0.0) {
	      SSL_ComplexAssign(-ra, -sa, &c1);
	      SSL_ComplexAssign(w, q, &c2);
	      SSL_ComplexDiv(c1, c2, &c3);
	      h[i-1][na-1]=c3.re;
	      h[i-1][en-1]=c3.im;
	    } /* if */
	    else {
	      x=h[i-1][i];
	      y=h[i][i-1];
	      vr=pow(wr[i-1]-p, 2.0)+wi[i-1]*wi[i-1]-q*q;
	      vi=(wr[i-1]-p)*2.0*q;
	      if ((vr==0.0) && (vi==0.0)) 
		vr=macheps*norm*(fabs(w)+fabs(q)+fabs(x)+fabs(y)+fabs(z));
	      SSL_ComplexAssign(x*r-z*ra+q*sa, x*s-z*sa-q*ra, &c1);
	      SSL_ComplexAssign(vr, vi, &c2);
	      SSL_ComplexDiv(c1, c2, &c3);
	      h[i-1][na-1]=c3.re;
	      h[i-1][en-1]=c3.im;
	      if (fabs(x)>(fabs(z)+fabs(q))) {
		h[i][na-1]=(-ra-w*h[i-1][na-1]+q*h[i-1][en-1])/x;
		h[i][en-1]=(-sa-w*h[i-1][en-1]-q*h[i-1][na-1])/x;
	      } /* if */
	      else {
		SSL_ComplexAssign(-r-y*h[i-1][na-1], -s-y*h[i-1][en-1],
			      &c1);
		SSL_ComplexAssign(z, q, &c2);
		SSL_ComplexDiv(c1, c2, &c3);
		h[i][na-1]=c3.re;
		h[i][en-1]=c3.im;
	      } /* else */
	    } /* else */
	  } /* else */
	} /* for i */
      } /* if */
  } /* for en */

  for(i=1; i<=(low-1); i++)
    for(j=(i+1); j<=n; j++)
      vecs[i-1][j-1]=h[i-1][j-1];
  for(i=(upp+1); i<=n; i++)
    for(j=(i+1); j<=n; j++) 
      vecs[i-1][j-1]=h[i-1][j-1];
  
  for(j=n; j>=low; j--) {
    if (j<=upp)
      m=j;
    else
      m=upp;
    l=j-1;
    if (wi[j-1]<0.0) {
      for(i=low; i<=upp; i++) {
	y=z=0.0;
	for(k=low; k<=m; k++) {
	  y+=vecs[i-1][k-1]*h[k-1][l-1];
	  z+=vecs[i-1][k-1]*h[k-1][j-1];
	} /* for k */
	vecs[i-1][l-1]=y;
	vecs[i-1][j-1]=z;
      } /* for i */
    } /* if */
    else
      if (wi[j-1]==0.0)
	for(i=low; i<=upp; i++) {
	  z=0.0;
	  for(k=low; k<=m; k++)
	    z+=vecs[i-1][k-1]*h[k-1][j-1];
	  vecs[i-1][j-1]=z;
	} /* for i */
  } /* for j */
  
 L270:
  dummy=0;
} /* hqr2 */
	  
void Eigen(int n, int ref, double **AJAC, int maxit, double eps, 
	   int fixedref, SSL_Complex *values, SSL_Complex **vectors) {

  double  *wr, *wi, *bald, **T, **A;
  int     i, j, ballow, balhi, max, block;
  int     *intout;
  int     fail;

  intout=IntVectorAlloc(n);
  wr=VectorAlloc(n);
  wi=VectorAlloc(n);
  bald=VectorAlloc(n);
  T=MatrixAlloc(n);
  A=MatrixAlloc(n);

  for(i=1; i<=n; i++)
    for(j=1; j<=n; j++)
      A[i-1][j-1]=AJAC[i-1][j-1];

  Balance(n, 10, A, &ballow, &balhi, bald);
  Elmhes(n, ballow, balhi, A, intout);
  Elmtrans(n, ballow, balhi, A, intout, T);

  hqr2(n, ballow, balhi, maxit, eps, A, T, wr, wi, intout, &fail);
  if (fail==1) 
    (void) fprintf(stderr, "Failure in hqr2 function. Do not trust the given eigenvectors and -values\n");
  /*
  tmxx=0;
  for(i=1; i<=n; i++)
    if (abs(intout[i-1])>tmxx)
      tmxx=(int)ceil(abs(intout[i-1]));
  */
  for(i=1; i<=n; i++)
    for(j=1; j<=n; j++)
      A[i-1][j-1]=0.0;
  i=1;
  do {
    if (wi[i-1]!=0.0) {
      A[i-1][i-1]=wr[i-1];
      A[i][i]=wr[i-1];
      A[i-1][i]=wi[i-1];
      A[i][i-1]=wi[i];
      i+=2;
    } /* if */
    else {
      A[i-1][i-1]=wr[i-1];
      i++;
    } /* else */
  } while (i<n);
  if (i==n)
    A[i-1][i-1]=wr[i-1];
 
  Swap(n, A, T, eps);
  BalBak(n, ballow, balhi, n, T, bald);
  NormalizingMatrix(n, A, fixedref, &ref, T, eps);

  /* store eigenvectors and eigenvalues nicely */
  i=1;              /* eigenvalues */
  do {
    BlockCheck(A, n, i, &block, eps);
    if (block==1) {
      SSL_ComplexAssign(A[i-1][i-1], A[i-1][i], &values[i-1]);
      SSL_ComplexAssign(A[i][i], A[i][i-1], &values[i]);
      i+=2;
    } else {
      SSL_ComplexAssign(A[i-1][i-1], 0.0, &values[i-1]);
      i++;
    } /* if else */
  } while (i!=(n+1));
  i=1;               /* eigenvectors */
  do {
    BlockCheck(A, n, i, &block, eps);
    if (block==1) {
      for(j=1; j<=n; j++)
	SSL_ComplexAssign(T[j-1][i-1], T[j-1][i], &vectors[i-1][j-1]);
      for(j=1; j<=n; j++)
	SSL_ComplexAssign(T[j-1][i-1], -T[j-1][i], &vectors[i][j-1]);
      i+=2;
    } else {
      for(j=1; j<=n; j++)
	SSL_ComplexAssign(T[j-1][i-1], 0.0, &vectors[i-1][j-1]);
      i++;
    } /* if else */
  } while (i!=(n+1));

  VectorFree(n, wi);
  VectorFree(n, wr);
  VectorFree(n, bald);
  IntVectorFree(n, intout);
  MatrixFree(n, A);
  MatrixFree(n, T);
} /* Eigen */