// Eigen2.java
// The procedure eigen computes both the eigenvalues and eigenvectors
// of arbitrary n by n complex matrix.
// This is a Java version of an Ada version of a NAG Fortran
// library subroutine.
// Some Fortran labels have been preserved for traceability

package myjava;
import myjava.*;

public strictfp class Eigen2
{
  public static void eigen(Complex A[][], Complex lambda[],
                           Complex vec[][], boolean fail[])
  {
    System.out.println("Eigene.eigen(A, lambda, vec, fail)");
    // driver for computing eigenvalues and eigenvectors
    if(A==null || lambda==null || vec==null)
    {
      System.out.println("Error in Eigen2.eigen,"+
                         " null or inconsistent array sizes.");
      return;
    }
    int n = A.length;
    if(A[0].length!=n || vec.length!=n || vec[0].length!=n || lambda.length!=n)
    {
      System.out.println("Error in Eigen.eigen,"+
                         " inconsistent array sizes.");
      return;
    }
    fail[0] = false;

    // special cases
    if(n<1) {System.out.println("zero size matrix"); return;}
    int rowcol[] = new int[n];
    Complex B[][] = new Complex[n][n];
    ComplexMatrix.copy(A, B);

    if(n==1)
    {
      lambda[0] = B[0][0];
      vec[0][0] = new Complex(1.0, 0.0);
      return;
    }
    if(n==2)
    {
      twobytwo(B, lambda, vec);
      return;
    }
    System.out.println("calling cxhess");
    cxhess(B, rowcol);
    for(int i=0; i<n; i++) lambda[i] = new Complex(-999.0, -999.0);
    System.out.println("calling cxeig2c");
    cxeig2c(B, lambda, vec, rowcol, fail);
  } // end eigen

  private static void twobytwo(Complex A[][], Complex lambda[],
                               Complex vec[][])
  {
    Complex b, c, rad, l1, l2;
    Complex Z[] = new Complex[2];
    double t;

    b = A[0][0].add(A[1][1]); // negative
    c = A[0][0].multiply(A[1][1]);
    c = c.subtract(A[0][1].multiply(A[1][0]));
    rad = (b.multiply(b)).subtract(c.multiply(4.0)); // a==1
    rad = rad.sqrt();
    l1 = (b.add(rad)).divide(2.0);
    l2 = (b.subtract(rad)).divide(2.0);
    lambda[0] = l1;
    lambda[1] = l2;
    // eigenvectors in columns
    Z[0] = A[0][1].negate();
    Z[1] = A[0][0].subtract(l1);
    t = ComplexMatrix.norm2(Z);
    vec[0][0] = Z[0].divide(t);
    vec[1][0] = Z[1].divide(t);
    Z[0] = A[1][1].subtract(l2);
    Z[1] = A[1][0].negate();
    t = ComplexMatrix.norm2(Z);
    vec[0][1] = Z[0].divide(t);
    vec[1][1] = Z[1].divide(t);
  }

  private static double sumabs(Complex Z)
  {
    return Math.abs(Z.real())+Math.abs(Z.imaginary());
  } // end sumabs

  private static void cxhess(Complex A[][], int rowcol[])
  {
    int i, k, t;
    Complex x;
    Complex y;

    System.out.println("cxhess(A, rowcol)");
    int n = A.length; // checked before call

    for(int j=0; j<n; j++) rowcol[j] = j;
    k = 0;
    for(int m=k+1; m<n-1; m++) // main reduction loop
    {
      i= m;
      x = new Complex(0.0, 0.0);
      for(int j=m; j<n; j++)
      {
	if(sumabs(A[j][m-1]) > sumabs(x))
        {
          x = A[j][m-1];
          i = j;
        }
      }
      if(i != m)
      {
        //  rowcol row column interchange of H
	t = rowcol[m];
        rowcol[m] = rowcol[i];
        rowcol[i] = t;
	for(int j=m-1; j<n; j++)
        {
          y = A[i][j];
          A[i][j] = A[m][j];
          A[m][j] = y;
        }
        for(int j=0; j<n; j++) //for J in 1..N loop
        {
          y = A[j][i];
          A[j][i] = A[j][m];
          A[j][m] = y;
        }
      }
      if(sumabs(x) != 0.0)
      {
	for(int ii=m+1; ii<n; ii++)
        {
          y = A[ii][m-1];
          if(sumabs(y) > 0.0)
	  {
            y = y.divide(x);
            A[ii][m-1] = y;
            for(int j=m; j<n; j++)
            {
              A[ii][j] = A[ii][j].subtract(y.multiply(A[m][j]));
	    }
            for(int j=0; j<n; j++)
            {
              A[j][m] = A[j][m].add(y.multiply(A[j][ii]));
            }
          } // end if
          A[ii][m-1] = new Complex(0.0, 0.0); // just cleanup
        }
      } // end if
    } // end main reduction loop
    System.out.println("result of cxhess=");
    ComplexMatrix.print(A);
    System.out.println("based on rowcol interchanges=");
    Matrix.print(rowcol);
    System.out.println(" ");
  } // end cxhess

  private static void cxeig2c(Complex A[][], Complex lambda[],
                              Complex vec[][],
                              int rowcol[], boolean fail[])
  {
    int j, k, m, mm, low, its, itn, ien;
    double anorm = 0.0;
    double ahr, aahr, acc, xr, xi, yr, yi, zr;
    Complex accnorm;
    Complex x, y, z, yy, T, S;
    int n = A.length; // checked in driver

    low = 0;
    acc = Math.pow(2.0, -23);
    System.out.println("acc="+acc+" = 2^-23");
    T = new Complex(0.0, 0.0);
    itn = 30 * n; // heuristic on maximum iterations
    ComplexMatrix.identity(vec); // initialize to identity Matrix

    // starting from Hessenberg reduction
    for(int ii=n-2; ii>0; ii--) //for i in reverse A'FIRST+1..A'LAST-1 loop
    {
      j = rowcol[ii];
      for(k=ii+1; k<n; k++) //for K in i+1..A'LAST loop
      {
        vec[k][ii] = A[k][ii-1];
      }
      if(ii != j)
      {
	for(k=ii; k<n; k++) //for k in i..A'LAST loop
        {
          vec[ii][k] = vec[j][k];
          vec[j][k] = new Complex(0.0, 0.0);
        }
        vec[j][ii] = new Complex(1.0, 0.0);
      }
    }
    ien = n-1; // used as subscript, loop test <=ien
 
    // ien is decremented
    while(low <= ien) // 260
    {
      System.out.println("260 low="+low+", ien="+ien);
      its = 0;
      // look for small single subdiagonal element
      L280: while(true)  //  280
      { 
        System.out.println("in 280");
        k = low;
        // for kk in reverse low+1..ien loop  // 300
        for(int kk=ien; kk>low; kk--) // 300
	{
          System.out.println("300 kk="+kk);

          ahr = sumabs(A[kk][kk-1]);
          aahr = acc * (sumabs(A[kk-1][kk-1]) + sumabs(A[kk][kk]));
          if(ahr <= aahr)
	  {
            k = kk;
            break;
          }
        }  //  300 
        System.out.println("exiting 300 with k="+k);
        if(k == ien) break L280;  //exit L280 when k = ien;  // 780
        if(itn <= 0)
	{
          fail[0] = true;
          return;
        }

        // compute shift
        if(its == 10 || its == 20)
	{
          S = new Complex(Math.abs(A[ien][ien-1].real()) +
                          Math.abs(A[ien-1][ien-2].real()) ,
                          Math.abs(A[ien][ien-1].imaginary()) +
                          Math.abs(A[ien-1][ien-2].imaginary()));
	}
        else
	{
          S = A[ien][ien];
          x = A[ien-1][ien].multiply(A[ien][ien-1]);
          if(sumabs(x) > 0.0)
	  {
            y = (A[ien-1][ien-1].subtract(S)).divide(new Complex(2.0, 0.0));
            z = ((y.multiply(y)).add(x)).sqrt();
            if(y.real() * z.real() + y.imaginary() * z.imaginary() < 0.0)
	    {
              z = z.negate();
            }
            yy = y.add(z);
            S = S.subtract(x.divide(yy));
          } // end if;
        } // end if;  //  400 
        for(int i=low; i<=ien; i++) // for i in low..ien loop  // 420
	{
          A[i][i] = A[i][i].subtract(S);
        } // end loop;  //  420 
        T = T.add(S);
        its = its + 1;
        itn = itn - 1;
        j = k + 1;

        // look for two consecutive small sub-diagonal elements
        xr = sumabs(A[ien-1][ien-1]);
        yr = sumabs(A[ien][ien-1]);
        zr = sumabs(A[ien][ien]);
        m = k;
        for(mm=ien-1; mm>=j; mm--) // for mm in reverse j..ien-1 loop  // 460
        {
	  System.out.println("460 mm="+mm+", m="+m+", j="+j);
          yi = yr;
          yr = sumabs(A[mm][mm-1]);
          xi = zr;
          zr = xr;
          xr = sumabs(A[mm-1][mm-1]);
          if(yr <= (acc * zr/yi *( zr + xr + xi )))
	  {
            m = mm;
            break;
          }
        } //end loop;  //  460 

        // triangular decomposition  A = L*R
        for(int i=m+1; i<=ien; i++) //for i in m+1..ien loop  // 620
        {
	  System.out.println("620 mm="+mm+", m="+m+", i="+i);
          x = A[i-1][i-1];
          y = A[i][i-1];
          if(sumabs(x) >= sumabs(y))
	  {
            z = y.divide(x);
            lambda[i] = new Complex(-1.0, 0.0 );
          }
          else
	  {
            // interchange rows of A
	    for(int jj=i-1; jj<n; jj++) // for j in i-1..n loop  // 540
            {
              z = A[i-1][jj];
              A[i-1][jj] = A[i][jj];
              A[i][jj] = z;
            } // end loop;  //  540
            z = x.divide(y);
            lambda[i] = new Complex(1.0, 0.0);
          } // end if;
          A[i][i-1] = z;
          for(int jj=i; jj<n; jj++) // for j in i .. N loop  // 600
          {
            A[i][jj] = A[i][jj].subtract(z.multiply(A[i-1][jj]));
          } // end loop;  //  600
        } // end loop;  //  620 

        // composition R*L = H
        for(int jj=m+1; jj<=ien; jj++) // for j in m+1..ien loop  // 760
        {
          x = A[jj][jj-1];
          A[jj][jj-1] = new Complex(0.0, 0.0);

          // interchange columns of A and vec if necessary
          if(lambda[jj].real() > 0.0)
          {
            for(int i=low; i<=jj; i++) // for i in low .. j loop  // 660
            {
              z = A[i][jj-1];
              A[i][jj-1] = A[i][jj];
              A[i][jj] = z;
            } // end loop;  //  660
            for(int i=low; i<n; i++) // for i in low .. N loop  // 680
            {
              z = vec[i][jj-1];
              vec[i][jj-1] = vec[i][jj];
              vec[i][jj] = z;
            } // end loop;  //  680
          } // end if

          // end interchange columns
          for(int i=low; i<=jj; i++) // for i in low..j loop  // 720
          {
            A[i][jj-1] = A[i][jj-1].add(x.multiply(A[i][jj]));
          }  //  720
          for(int i=low; i<n; i++) // for i in low..N loop  // 740
          {
            vec[i][jj-1] = vec[i][jj-1].add(x.multiply(vec[i][jj]));
          }  //  740

          // end accumulate transformations
        }  //  760 
      }  // 280
 
      // a root found
      lambda[ien] = A[ien][ien].add(T);
      ien = ien - 1;
    } // end loop;  // 260 while

    // all roots found
    for(int i=0; i<n; i++) // for i in A'RANGE loop
    {
      anorm = anorm + sumabs(lambda[i]);
      for(int jj=i+1; jj<n; jj++) // for j in i + 1 .. A'LAST loop
      {
        anorm = anorm + sumabs(A[i][jj]);
      }
    }
    accnorm = new Complex(anorm * Math.pow(2.0,-23), 0.0);
    if(anorm == 0.0 || n < 2)
    {
      return;  // done
    }

    // back substitute to set up vec of upper triangular form
    for(ien=n-1; ien>low; ien--) // for ien in reverse low+1..N loop
    {
      x = lambda[ien];
      for(int i=ien-1; i>=low; i--) // for i in reverse low .. ien - 1 loop
      {
        z = A[i][ien];
        for(int jj=i+1; jj<ien; jj++) // for j in i+1..ien-1 loop
        {
          z = z.add(A[i][jj].multiply(A[jj][ien]));
        }
        y = x.subtract(lambda[i]);
        if(sumabs(y) == 0.0)
        {
          y = accnorm;
        }
        A[i][ien] = z.divide(y);
      }
    }

    // multiply by transformation Matrix to give vec of original full Matrix
    for(int jj=n-1; jj>=0; jj--) // for j in reverse A'RANGE loop
    {
      for(int i=0; i<n; i++) // for i in A'RANGE loop
      {
        z = vec[i][jj];
        for(k=0; k<jj; k++) // for k in A'first..j-1 loop
	{
          z = z.add(vec[i][k].multiply(A[k][jj]));
        }
        vec[i][jj] = z;
      }
    }
  } // end cxeig2c
} // end class Eigen2