// TestEigen2  just developmental test

package myjava;
import myjava.*;

class TestEigen2
{
  public static void main(String args[])
  {
    boolean fail[] = new boolean[1];
    Complex A1[][] = new Complex[1][1];
    Complex VEC1[][] = new Complex[1][1];
    Complex VAL1[] = new Complex[1];

    Complex A2[][] = new Complex[2][2];
    Complex B2[][] = new Complex[2][2];
    Complex X2[] = new Complex[2];
    Complex Y2[] = new Complex[2];
    Complex VEC2[][] = new Complex[2][2];
    Complex VAL2[] = new Complex[2];

    Complex A3[][] = new Complex[3][3];
    Complex VEC3[][] = new Complex[3][3];
    Complex VAL3[] = new Complex[3];

    Complex A4[][] = new Complex[4][4];
    Complex VEC4[][] = new Complex[4][4];
    Complex VAL4[] = new Complex[4];

    A1[0][0] = new Complex(2.0,3.0);
    System.out.println("initial A1");
    ComplexMatrix.print(A1);
    System.out.println("eigen2(A1, VAL1, VEC1, fail");
    Eigen2.eigen(A1, VAL1, VEC1, fail);
    System.out.println("fail="+fail[0]+",  VEC1=");
    ComplexMatrix.print(VEC1);
    System.out.println("VAL1=");
    ComplexMatrix.print(VAL1);
    System.out.println("eigenCheck(A1, VEC1, VAL1)");
    ComplexMatrix.eigenCheck(A1, VEC1, VAL1);
    System.out.println();

    A2[0][0] = new Complex(1.0,1.0);
    A2[0][1] = new Complex(2.0,-3.0);
    A2[1][0] = new Complex(5.0,-1.0);
    A2[1][1] = new Complex(2.0,2.0);
    System.out.println("initial A2");
    ComplexMatrix.print(A2);
    System.out.println("eigen2(A2, VAL2, VEC2, fail");
    Eigen2.eigen(A2, VAL2, VEC2, fail);
    System.out.println("fail="+fail[0]+",  VEC2=");
    ComplexMatrix.print(VEC2);
    System.out.println("VAL2=");
    ComplexMatrix.print(VAL2);
    B2[0][0] = A2[0][0].subtract(VAL2[0]);
    B2[1][0] = A2[1][0];
    B2[0][1] = A2[0][1];
    B2[1][1] = A2[1][1].subtract(VAL2[0]);
    System.out.println("solve |A-lambda I||v|=|0|");
    ComplexMatrix.print(B2);
    System.out.println("det |A-lambda I|="+
                    ComplexMatrix.determinant(B2));
    Y2[0] = new Complex(0.0, 0.0);
    Y2[1] = new Complex(0.0, 0.0);
    ComplexMatrix.solve(B2, Y2, X2);
    System.out.println("first eigenvector=");
    ComplexMatrix.print(X2);
    B2[0][0] = A2[0][0].subtract(VAL2[1]);
    B2[1][0] = A2[1][0];
    B2[0][1] = A2[0][1];
    B2[1][1] = A2[1][1].subtract(VAL2[1]);
    ComplexMatrix.print(B2);
    System.out.println("det |A-lambda I|="+
                    ComplexMatrix.determinant(B2));
    Y2[0] = new Complex(0.0, 0.0);
    Y2[1] = new Complex(0.0, 0.0);
    ComplexMatrix.solve(B2, Y2, X2);
    System.out.println("second eigenvector=");
    ComplexMatrix.print(X2);

    System.out.println("eigenCheck(A2, VEC2, VAL2)");
    ComplexMatrix.eigenCheck(A2, VEC2, VAL2);
    System.out.println();

    A3[0][0] = new Complex(1.0,1.0);
    A3[0][1] = new Complex(2.0,-3.0);
    A3[0][2] = new Complex(5.0,-1.0);
    A3[1][0] = new Complex(2.0,-3.0);
    A3[1][1] = new Complex(1.0,1.0);
    A3[1][2] = new Complex(2.0,-5.0);
    A3[2][0] = new Complex(5.0,-1.0);
    A3[2][1] = new Complex(2.0,-5.0);
    A3[2][2] = new Complex(4.0,-4.0);
    System.out.println("initial A3");
    ComplexMatrix.print(A3);
    System.out.println("eigen2(A3, VAL3, VEC3, fail");
    Eigen2.eigen(A3, VAL3, VEC3, fail);
    System.out.println("fail="+fail[0]+",  VEC3=");
    ComplexMatrix.print(VEC3);
    System.out.println("VAL3=");
    ComplexMatrix.print(VAL3);
    System.out.println("eigenCheck(A3, VEC3, VAL3)");
    ComplexMatrix.eigenCheck(A3, VEC3, VAL3);
    System.out.println("trying eigvec");
    eigvec(A3, VAL3, VEC3);
    System.out.println("VEC3=");
    ComplexMatrix.print(VEC3);
    System.out.println("eigenCheck(A3, VEC3, VAL3)");
    ComplexMatrix.eigenCheck(A3, VEC3, VAL3);
    System.out.println();

    A4[0][0] = new Complex(1.0,4.0);
    A4[1][0] = new Complex(2.0,4.0);
    A4[2][0] = new Complex(3.0,4.0);
    A4[3][0] = new Complex(4.0,4.0);
    A4[0][1] = new Complex(2.0,4.0);
    A4[1][1] = new Complex(2.0,3.0);
    A4[2][1] = new Complex(3.0,3.0);
    A4[3][1] = new Complex(4.0,3.0);
    A4[0][2] = new Complex(3.0,4.0);
    A4[1][2] = new Complex(3.0,3.0);
    A4[2][2] = new Complex(3.0,2.0);
    A4[3][2] = new Complex(4.0,2.0);
    A4[0][3] = new Complex(4.0,4.0);
    A4[1][3] = new Complex(4.0,3.0);
    A4[2][3] = new Complex(4.0,2.0);
    A4[3][3] = new Complex(4.0,1.0);
    System.out.println("initial A4");
    ComplexMatrix.print(A4);
    System.out.println("eigen2(A4, VAL4, VEC4, fail");
    Eigen2.eigen(A4, VAL4, VEC4, fail);
    System.out.println("fail="+fail[0]+",  VEC4=");
    ComplexMatrix.print(VEC4);
    System.out.println("VAL4=");
    ComplexMatrix.print(VAL4);
    System.out.println("eigenCheck(A4, VEC4, VAL4)");
    ComplexMatrix.eigenCheck(A4, VEC4, VAL4);
    System.out.println();
    System.out.println("trying eigvec");
    eigvec(A4, VAL4, VEC4);
    System.out.println("VEC4=");
    ComplexMatrix.print(VEC4);
    System.out.println("eigenCheck(A4, VEC4, VAL4)");
    ComplexMatrix.eigenCheck(A4, VEC4, VAL4);
    System.out.println();
  }

  public static void eigvec(final Complex A[][], final Complex Y[],
                            Complex X[][])
  {
    // Y[i] are eigenvalues of A, X[i][] is eigenvector
    int n=A.length;
    int m=n+1;
    Complex B[][]=new Complex[n][m];   // working matrix
    int row[]=new int[n];              // row interchange indicies
    int hold , I_pivot;                // pivot indicies
    Complex pivot;                     // pivot element value
    double abs_pivot;
    Complex T[] = new Complex[n];
    double r;
    Complex z;

    if(A[0].length!=n || Y.length!=n || X.length!=n)
    {
      System.out.println("Error in simeq,"+
                         " inconsistent array sizes.");
    }
    // build working data structure
    for(int v=0; v<n; v++) // loop through eigenvalues
    {
      System.out.println("working eigenvalue v="+v);
      for(int i=0; i<n; i++)
      {
        for(int j=0; j<n; j++)
        {
          B[i][j] = A[i][j];
          if(i==j) B[i][i] = B[i][i].subtract(Y[v]);
        }
        B[i][n] = new Complex(0.0, 0.0);
      }
      System.out.println("initial B for eigenvalue Y["+v+"]="+Y[v]);
      ComplexMatrix.print(B);
      // set up row interchange vectors
      for(int k=0; k<n; k++)
      {
        row[k] = k;
      }
      //  begin main reduction loop
      for(int k=0; k<n; k++)
      {
        // find largest element for pivot
        pivot = B[row[k]][k];
        abs_pivot = pivot.abs();
        I_pivot = k;
        for(int i=k+1; i<n; i++)
        {
          if(B[row[i]][k].abs() > abs_pivot)
          {
            I_pivot = i;
            pivot = B[row[i]][k];
            abs_pivot = pivot.abs();
          }
        }
        // have pivot, interchange row indicies
        hold = row[k];
        row[k] = row[I_pivot];
        row[I_pivot] = hold;
	System.out.println("k="+k+", row[k]="+row[k]);
        // check for near singular
        if(abs_pivot < 1.0E-6)
        {
          for(int j=k; j<n+1; j++)
          {
            B[row[k]][j] = new Complex(0.0, 0.0);
          }
          System.out.println("redundant row, break for vector computation");
          break;
        } 
        else
        {
          // reduce about pivot
	    for(int j=k+1; j<n; j++)
          {
            B[row[k]][j] = B[row[k]][j].divide(B[row[k]][k]);
          }
	  B[row[k]][k] = new Complex(1.0, 0.0);
          //  inner reduction loop
          for(int i=k+1; i<n; i++)
          {
            for(int j=k+1; j<n; j++)
            {
              B[row[i]][j] = B[row[i]][j].subtract(
                             B[row[i]][k].multiply(B[row[k]][j]));
            }
            B[row[i]][k] = new Complex(0.0, 0.0);
          }
        }
        //  finished inner reduction
      }
      //  end main reduction loop
      //  backsubstitute  X  for return, unscrambling rows
      System.out.println("backsubstitute B for Y["+v+"]="+Y[v]);
      ComplexMatrix.print(B);
      // check if last row, simple back substitute
      // else  1, 0, rest  then 0, 1, rest etc
      B[row[n-1]][n] = new Complex(1.0, 0.0);
      X[n-1][v] = B[row[n-1]][n];
      T[n-1] = B[row[n-1]][n];
      for(int i=n-2; i>=0; i--)
      {
	System.out.println("back row i="+i);
	z = new Complex(0.0, 0.0);
	for(int jj=i+1; jj<n; jj++)
        {
	  System.out.println("back col jj="+jj);
	  z = z.add(B[row[i]][jj].multiply(B[row[jj]][n]));
        }
        B[row[i]][n] = (z.divide(B[row[i]][i])).negate();
        X[i][v] = B[row[i]][n];
        System.out.println("B["+row[i]+"]["+n+"] ="+B[row[i]][n]); 
        T[i] = B[row[i]][n];
      }
      r = ComplexMatrix.norm2(T);
      System.out.println("norm2 of X ="+r);
      for(int i=0; i<n; i++)
      {
        X[i][v] = X[i][v].divide(r);
        System.out.println("X["+i+"]["+v+"] ="+X[i][v]); 
      }
    }
  } // end eigvec
} // end file TestEigen2.java