/*				 16.FEB.2000			   tqli_mat.cc
						 from "Numerical Recipes in C"

  On input d[0...n-1] contains the diagonal elements of the tridiagonal
  matrix, and e[0...n-2] the subdiagonal elements. If the eigenvectors
  of a tridiagonal matrix are desired, z[0 .. n-1][0 .. n-1] must be the 
  identity matrix, otherwise z is the output transformation matrix of tred2.
  On output d contains the eigenvalues, e is destroyed and the columns
  of z contain the eigenvectors.
  This function returns 1 if succesful, 0 otherwise.

*/

#include <stdio.h>
#include <math.h>

#include "corso.h"

#define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))

int tqli(double *d,double *e,double **z,int n)
{
  int iter,i,k,l,m;
  double b,c,dd,f,g,p,r,s;

  for (l=0;l<n;l++)
    { iter=0;
      do
	{ for (m=l;m<n-1;m++)
	    { dd=fabs(d[m])+fabs(d[m+1]);
	      if (fabs(e[m])+dd == dd) break;
	    }
	  if (m!=l)
	    { if (iter++ == 30)
		{ printf("Too many iterations in tqli()!!!\n");
		  return 0;
		}
	      g=(d[l+1]-d[l])/(2.0*e[l]);
	      r=sqrt((g*g)+1.0);
	      g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
	      s=c=1.0;
	      p=0.0;
	      for (i=m-1;i>=l;i--)
		{ f=s*e[i];b=c*e[i];
		  if (fabs(f) >= fabs(g))
		    { c=g/f;r=sqrt((c*c)+1.0);
		      e[i+1]=f*r;
		      c *= (s=1.0/r);
		    }
		  else
		    { s=f/g;r=sqrt((s*s)+1.0);
		      e[i+1]=g*r;
		      s *= (c=1.0/r);
		    }
		  g=d[i+1]-p;
		  r=(d[i]-g)*s+2.0*c*b;
		  p=s*r;
		  d[i+1]=g+p;
		  g=c*r-b;
		  for (k=0;k<n;k++)
		    { f=z[k][i+1];
		      z[k][i+1]=s*z[k][i]+c*f;
		      z[k][i]=c*z[k][i]-s*f;
		    }
		}
	      d[l]=d[l]-p;
	      e[l]=g;
	      e[m]=0.0;
	    }
	} while (m!=l);
    }
  return 1;
}