File Coverage

lib/PDL/MatrixOps/eigen.c

Criterion	Covered	Total	%
statement	493	633	77.8
branch	250	328	76.2
condition			n/a
subroutine			n/a
pod			n/a
total	743	961	77.3

/*
 * eigen.c - calculation of eigen values and vectors.
 *
 * (C) Copyright 2001 by NetGroup A/S. All rights reserved.
 *
 * $Log$
 * 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
#include "matrix.h"
#include
#include

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

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

  if (i>=n-1)
    *block=0;
  else {
    if ((fabs(A[i][i+1]-A[i+1][i])>epsx) &&
	(fabs(A[i][i]-A[i+1][i+1])<=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=0;
  do {
    BlockCheck(A, n, i, &block, eps);
    if (block==1) {
      fprintf(outfile, "\t\t\t\t%e\t\t%e\n", A[i][i], A[i][i+1]);
      fprintf(outfile, "\t\t\t\t%e\t\t%e\n", A[i+1][i+1], A[i+1][i]);
      i+=2;
    } else {
      fprintf(outfile, "\t\t\t\t%e\t\t%e\n", A[i][i], 0.0);
      i++;
    } /* if else */
  } while (i
  fprintf(outfile, "\nEigenvectors:\t\t\tRe\t\t\tIm\n");
  i=0;
  do {
    BlockCheck(A, n, i, &block, eps);
    if (block==1) {
      for(j=0; j
	fprintf(outfile, "\t\t\t\t%e\t\t%e\n", B[j][i], B[j][i+1]);
      fprintf(outfile, "\n");
      for(j=0; j
	fprintf(outfile, "\t\t\t\t%e\t\t%e\n", B[j][i], -B[j][i+1]);
      fprintf(outfile, "\n");
      i+=2;
    } else {
      for(j=0; j
	fprintf(outfile, "\t\t\t\t%e\t\t%e\n", B[j][i], 0.0);
      fprintf(outfile, "\n");
      i++;
    } /* if else */
  } while (i
} /* PrintEigen */


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

  int      j, col=0, block;

  do {
    double sumsq = 0;
    BlockCheck(A, n, col, &block, eps);
    for(j=0; j
      sumsq += (V[j][col] * V[j][col]) + (block==1 ? (V[j][col+1] * V[j][col+1]) : 0);
    double norm = sqrt(sumsq);
    if (norm == 0.0) continue;
    if (block==1) {
      for(j=0; j
	complex double c2 = V[j][col] + I * V[j][col+1], c3 = c2 / norm;
	V[j][col]=creal(c3);
	V[j][col+1]=cimag(c3);
      } /* for j */
      col+=2;
    } else {
      for(j=0; j
	V[j][col] /= norm;
      col++;
    }
  } while (col
} /* NormalizingMatrix */

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

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

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

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

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

  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=0;
  k=n-1;

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

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

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

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

    L140:
    if (c
      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]*=f;
      noconv=1;
      for(j=l; j
	a[i][j]*=g;
      for(j=0; j<=k; j++)
	a[j][i]*=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
    s=d[i];
    for(j=0; j
      z[i][j]*=s;
  } /* for i */
  for(i=low-1; i>=0; i--) {
    k=(int)floor(d[i]+0.5)-1;
    if (k!=i)
      for(j=0; j
	s=z[i][j];
	z[i][j]=z[k][j];
	z[k][j]=s;
      } /* for j */
  } /* for i */
  for(i=hi+1; i
    k=(int)floor(d[i]+0.5)-1;
    if (k!=i)
      for(j=0; j
	s=z[i][j];
	z[i][j]=z[k][j];
	z[k][j]=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;
  for(m=k+1; m
    i=m;
    x=0.0;
    for(j=m; j
      if (fabs(a[j][m-1])>fabs(x)) {
	x=a[j][m-1];
	i=j-1;
      } /* if */
    index[m]=i+1;
    if (i!=m) {
      for(j=m-1; j
        y=a[i][j];
        a[i][j]=a[m][j];
        a[m][j]=y;
      } /* for j */
      for(j=0; j
        y=a[j][i];
        a[j][i]=a[j][m];
        a[j][m]=y;
      } /* for j */
    } /* if */
    if (x!=0.0)
      for(i=(m+1); i
        y=a[i][m-1];
        if (y!=0.0) {
          a[i][m-1]=y/x;
          y/=x;
          for(j=m; j
            a[i][j]-=y*a[m][j];
          for(j=0; j
            a[j][m]+=y*a[j][i];
        } /* 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=0; i
    for(j=0; j
      v[i][j]=0.0;
    v[i][i]=1.0;
  } /* for i */
  for(i=(upp-1); i>=low+1; i--) {
    j=index[i]-1;
    for(k=(i+1); k
      v[k][i]=h[k][i-1];
    if (i!=j) {
      for(k=i; k
	v[i][k]=v[j][k];
	v[j][k]=0.0;
      } /* for k */
      v[j][i]=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;
  double  p = 0, q = 0, r = 0, s = 0, t, w, x, y, z = 0, ra, sa, vr, vi, norm;
  int     notlast;
  complex double c1, c2, c3;

  *fail=0;
  for(i=0; i
    wr[i]=h[i][i];
    wi[i]=0.0;
    cnt[i]=0;
  } /* for i */
  for(i=upp+1; i
    wr[i]=h[i][i];
    wi[i]=0.0;
    cnt[i]=0;
  } /* for i */
  en=upp;
  t=0.0;

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

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

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

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

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

  for(k=m; k
    if (k!=na)
      notlast=1;
    else
      notlast=0;
    if (k!=m) {
      p=h[k][k-1];
      q=h[k+1][k-1];
      if (notlast==1)
	r=h[k+2][k-1];
      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+1!=m+1)
      h[k][k-1]=-s*x;
    else
      if (l!=m+1)
	h[k][k-1]=-h[k][k-1];
    p+=s;
    x=p/s;
    y=q/s;
    z=r/s;
    q/=p;
    r/=p;

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

    for(i=0; i
      p=x*h[i][k]+y*h[i][k+1];
      if (notlast==1) {
	p+=z*h[i][k+2];
	h[i][k+2]-=p*r;
      } /* if */
      h[i][k+1]-=p*q;
      h[i][k]-=p;
    } /* for i */

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

  L233:;
  } /* for k */
  goto L220;

  L240:
  h[en][en]=x+t;
  wr[en]=h[en][en];
  wi[en]=0.0;
  cnt[en]=its;
  en=na;
  goto L210;

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

      /* h[en][j]=q*h[en][j]-p*z */
      h[en][j]*=q;
      h[en][j]-=p*z;
    } /* for j */
    for(i=0; i
      z=h[i][na];
      h[i][na]=q*z+p*h[i][en];

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

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

  L260:
  norm=0.0;
  k=0;
  for(i=0; i
    for(j=k; j
      norm+=fabs(h[i][j]);
    k=i;
  } /* for i */

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

  for(i=0; i
    for(j=i+1; j
      vecs[i][j]=h[i][j];
  for(i=upp+1; i
    for(j=i-1; j
      vecs[i][j]=h[i][j];

  for(j=n-1; j>=low; j--) {
    if (j<=upp)
      m=j+1;
    else
      m=upp+1;
    l=j-1;
    if (wi[j]<0.0) {
      for(i=low; i
	y=z=0.0;
	for(k=low; k
	  y+=vecs[i][k]*h[k][l];
	  z+=vecs[i][k]*h[k][j];
	} /* for k */
	vecs[i][l]=y;
	vecs[i][j]=z;
      } /* for i */
    } /* if */
    else
      if (wi[j]==0.0)
	for(i=low; i
	  z=0.0;
	  for(k=low; k
	    z+=vecs[i][k]*h[k][j];
	  vecs[i][j]=z;
	} /* for i */
  } /* for j */

 L270:;
} /* hqr2 */

char *Eigen(int n, double *AJAC, int maxit, double eps,
	   complex double *values, complex double *vectors) {

  double  *wr, *wi, *bald, **T, **A;
  int     i, j, ballow, balhi, 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=0; i
    for(j=0; j
      A[i][j]=AJAC[i*n + j];

  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)
    return "Failure in hqr2 function";
  for(i=0; i
    for(j=0; j
      A[i][j]=0.0;
  i=0;
  do {
    if (wi[i]!=0.0) {
      A[i][i]=wr[i];
      A[i+1][i+1]=wr[i];
      A[i][i+1]=wi[i];
      A[i+1][i]=wi[i+1];
      i+=2;
    } /* if */
    else {
      A[i][i]=wr[i];
      i++;
    } /* else */
  } while (i
  if (i==n-1)
    A[i][i]=wr[i];

  Swap(n, A, T, eps);
  BalBak(n, ballow, balhi, n, T, bald);
  NormalizingMatrix(n, A, T, eps);

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

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