import java.io.*;
import javax.swing.*;


class NA135
{
  //===========Eigen value calculations ==============

public static double[][] balance(double b[][])
{
// balance of a matrix for more accurate eigenvalue
// calculations
double radix=2.0;
double sqrdx=radix*radix;
double c,r,f,s,g;
int m,j,i,last;
int n=b.length;
last=0;
double a[][];
a=new double[n][n];
f=1;
s=1;
for(i=1;i<=n;i++)
  for(j=1;j<=n;j++)
    a[i-1][j-1]=b[i-1][j-1];
while(last==0)
  {
  last=1;
  for(i=1;i<=n;i++)
    {
    c=0;r=0;
    for(j=1;j<=n;j++)
      {
      if(j != i)
	{
	c+=Math.abs(a[j-1][i-1]);
	r+=Math.abs(a[i-1][j-1]);
	} //end of if(j!=..
      } //end of for(j=1...
    if(c != 0 &&  r != 0 )
      {
      g=r/radix;
      f=1.0;
      s=c+r;
      while(c<g)
	{
	f*=radix;
	c*=sqrdx;
	}
      g=r*radix;
      while(c>g)
	{
	f/=radix;
	c/=sqrdx;
	}
      } //end of if(c != 0 && ....
    if( (c+r)/f < 0.95*s )
      {
      last=0;
      g=1.0/f;
      for(j=1;j<=n;j++)  { a[i-1][j-1]*=g; }
      for(j=1;j<=n;j++)  { a[j-1][i-1]*=f; }
      } //end of if( ((c+r..
    }//end of for(i=1;i<=n....
  } //end of while last==0
return a;
}

public static double[][] Hessenberg(double b[][])
{
// Calculates the hessenberg matrix
// it is used in QR method to calculate eigenvalues
// of a matrix(symmetric or non-symmetric)
int m,j,i;
int n=b.length;
double a[][];
a=new double[n][n];
for(i=0;i<n;i++)
  for(j=0;j<n;j++)
    a[i][j]=b[i][j];
double x,y;
if(n>2)
  {
  for(m=2;m<=(n-1);m++)
     {
     x=0.0;
     i=m;
     for(j=m;j<=n;j++)
       {
       if(Math.abs(a[j-1][m-2]) > Math.abs(x))
	  {
	  x=a[j-1][m-2];
	  i=j;
	  } //end of if(Math.abs(..
       }//end of for(j=m,j<=n...
     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;
	 }//end of for(j=(m-1)..
       for(j=1;j<=n;j++)
	 {
	 y=a[j-1][i-1];
	 a[j-1][i-1]=a[j-1][m-1];
	 a[j-1][m-1]=y;
	 }//end of for(j=1;j<=n....
       } //end of if(i!=m)
     if(x != 0.0)
       {
       for(i=(m+1);i<=n;i++)
	  {
	  y=a[i-1][m-2];
	  if(y!=0.0)
	    {
	    y=y/x;
	    a[i-1][m-2]=y;
	    for(j=m;j<=n;j++)
	       {
	       a[i-1][j-1]-=y*a[m-1][j-1];
	       }
	    for(j=1;j<=n;j++)
	       {
	       a[j-1][m-1]+=y*a[j-1][i-1];
	       }
	    }//end of if(y!=0..
	  }//end of for(i=(m+1)...
       } //end of if(x != 0.0...
     }//end of for(m=2;m<=(n-1)..
  }//end of Hessenberg
for(i=1;i<=n;i++)
  for(j=1;j<=n;j++)
  {
  if(i>(j+1)) a[i-1][j-1]=0;
  }
return a;
}

public static double[][] QR(double b[][])
{
//calculates eigenvalues of a Hessenberg matrix
int n=b.length;
double rm[][]=new double[2][n];
double a[][]=new double[n+1][n+1];
double wr[]=new double[n+1];
double wi[]=new double[n+1];
int nn,m,l,k,j,its,i,mmin;
double z,y,x,w,v,u,t,s,r=0,q=0,p=0,anorm;
for(i=0;i<n;i++)
   for(j=0;j<n;j++)
      a[i+1][j+1]=b[i][j];
anorm=Math.abs(a[1][1]);
for (i=2;i<=n;i++)
   for (j=(i-1);j<=n;j++)
       anorm += Math.abs(a[i][j]);
nn=n;
t=0.0;
while (nn >= 1) {
   its=0;
   do {
      for (l=nn;l>=2;l--) {
      s=Math.abs(a[l-1][l-1])+Math.abs(a[l][l]);
      if (s == 0.0) s=anorm;
      if ((Math.abs(a[l][l-1]) + s) == s) break;
                          }
      x=a[nn][nn];
      if (l == nn) {
          wr[nn]=x+t;
	  wi[nn--]=0.0;
		   }
      else {
           y=a[nn-1][nn-1];
	   w=a[nn][nn-1]*a[nn-1][nn];
	   if (l == (nn-1)) {
	      p=0.5*(y-x);
	      q=p*p+w;
	      z=Math.sqrt(Math.abs(q));
	      x += t;
	      if (q >= 0.0) {
	         z=p+Matrix.SIGN(z,p);
		 wr[nn-1]=wr[nn]=x+z;
		 if (z!=0) wr[nn]=x-w/z;
		 wi[nn-1]=wi[nn]=0.0;
			    }
	      else {
	         wr[nn-1]=wr[nn]=x+p;
		 wi[nn-1]= -(wi[nn]=z);
		   }
	      nn -= 2;
			     }
	   else {
	      if (its == 30) System.out.println("Matrix.java, module QR, Too many iterations in hqr");
	      if (its == 10 || its == 20) {
	         t += x;
		 for (i=1;i<=nn;i++) a[i][i] -= x;
		 s=Math.abs(a[nn][nn-1])+Math.abs(a[nn-1][nn-2]);
		 y=x=0.75*s;
		 w = -0.4375*s*s;
					  }
	      ++its;
	      for (m=(nn-2);m>=l;m--) {
	      z=a[m][m];
	      r=x-z;
	      s=y-z;
	      p=(r*s-w)/a[m+1][m]+a[m][m+1];
	      q=a[m+1][m+1]-z-r-s;
	      r=a[m+2][m+1];
	      s=Math.abs(p)+Math.abs(q)+Math.abs(r);
	      p /= s;
	      q /= s;
	      r /= s;
	      if (m == l) break;
	      u=Math.abs(a[m][m-1])*(Math.abs(q)+Math.abs(r));
	      v=Math.abs(p)*(Math.abs(a[m-1][m-1])+
            Math.abs(z)+Math.abs(a[m+1][m+1]));
	      if ((u+v) == v) break;
					}
	      for (i=m+2;i<=nn;i++) {
	         a[i][i-2]=0.0;
	         if (i != (m+2)) a[i][i-3]=0.0;
		                    }
	      for (k=m;k<=nn-1;k++) {
	      if (k != m) {
	         p=a[k][k-1];
		 q=a[k+1][k-1];
		 r=0.0;
		 if (k != (nn-1)) r=a[k+2][k-1];
		 if ((x=Math.abs(p)+Math.abs(q)+Math.abs(r)) != 0.0) {
		    p /= x;
		    q /= x;
		    r /= x;
		                                                     }		                    }
		 if ((s=Matrix.SIGN(Math.sqrt(p*p+q*q+r*r),p)) != 0.0) {
		    if (k == m) {
		       if (l != m)
		       a[k][k-1] = -a[k][k-1];
				}
		    else
		       a[k][k-1] = -s*x;
		    p += s;
		    x=p/s;
		    y=q/s;
		    z=r/s;
		    q /= p;
		    r /= p;
		    for (j=k;j<=nn;j++) {
		    p=a[k][j]+q*a[k+1][j];
		    if (k != (nn-1)) {
		       p += r*a[k+2][j];
		       a[k+2][j] -= p*z;
		                     }
		    a[k+1][j] -= p*y;
		    a[k][j] -= p*x;
				       }
		    mmin = nn<k+3 ? nn : k+3;
		    for (i=l;i<=mmin;i++) {
		       p=x*a[i][k]+y*a[i][k+1];
		       if (k != (nn-1)) {
		          p += z*a[i][k+2];
			  a[i][k+2] -= p*r;
			                }
		       a[i][k+1] -= p*q;
		       a[i][k] -= p;
		  }
	       }
	    }
         }
      }
   } while (l < nn-1);
}
for(i=0;i<n;i++)
{
rm[0][i]=wr[i+1];
rm[1][i]=wi[i+1];
}
return rm;
} //end of QR

public static double[][] eigenValue(double b[][])
{
// this routine input a matrix (non symetric or symmetric)
// and calculate eigen values
// method balance can be used prior to this method to balance
// the input matrix
int n=b.length;
double d[][]=new double[2][n];
d=QR(Hessenberg(b));
return d;
}

public static complex[] eigenValueC(double b[][])
{
// this routine input a matrix (non symetric or symmetric)
// and calculate eigen values
// method balance can be used prior to this method to balance
// the input matrix
//output eigenvalues will be in a vector of complex form
int n=b.length;
double d[][]=new double[2][n];
d=QR(Hessenberg(b));
complex c[]=new complex[n];
for(int i=0;i<n;i++)
  {
  c[i]=new complex(d[0][i],d[1][i]);
  }
return c;
}
  public static void main(String args[]) throws IOException
  {
  //Eigen Value calculations using  method
  //double a[][]={{3,-1,0},{-1,2,-1},{0,-1,3}};
  //double a[][]={{1,-2,3},{-2,4,1},{3,1,2}};
  double a[][]={{3.556,-1.778,0},{-1.778,3.556,-1.778},{0,-1.778,3.556}};
  String s="Original matrix : ";
  s+=Matrix.toString(a);
  s+="QL solution matrix \n";
  s+=Matrix.toString(SCO13A.eigenQL(a));
  //Eigenvlue calculation
  complex c1[];
  int n=a.length;
  c1=new complex[n];
  c1=eigenValueC(a);
  s+="QR method eigenvalues \n";
  s+=Matrix.toStringT(c1);
  JOptionPane.showMessageDialog(null,s,
     "QR - QL eigenvalue method : ",JOptionPane.PLAIN_MESSAGE);
  }
}
