/*
 * Open Source Physics software is free software as described near the bottom of this code file.
 *
 * For additional information and documentation on Open Source Physics please see:
 * <http://www.opensourcephysics.org/>
 */

/**
    Adapted to OSP by Javier E. Hasbun and Wolfgang Christian, 2009.
  <P>
    @Copyright (c) 2009
    This software is to support the Open Source Physics library
    http://www.opensourcephysics.org under the terms of the GNU General Public
    License (GPL) as published by the Free Software Foundation.
 **/
package org.opensourcephysics.numerics;

/**
 * A Java Complex Eigenvalue Decomposition based on an Ada version of a NAG Fortran library subroutine.
 * Some Fortran labels have been preserved for traceability
 */
public class ComplexEigenvalueDecomposition implements java.io.Serializable {
  /**
   * Constructs the EigenvalueDecomposition.
   */
  public ComplexEigenvalueDecomposition() {}

  /**
   * Obtains the eigenvalues and eigenvectors of a complex matrix.
   *
   * Input
   * @parameter Complex double [][] matrix A containing the data
   * @parameter Complex double [] empty vector lambda
   * @parameter Complex double [][] empty eigenvector matrix vec
   * @parameter unset single element boolean [] fail
   *
   * Ouput:
   * matrix A is unmodified
   * vector lambda contains the complex eigenvalues
   * matrix vec contains the complex eigenvectors
   * parameter fails is set to true if the eigen fails to succeed
   * The procedure eigen computes both the eigenvalues and eigenvectors
   * of arbitrary n by n complex matrix.
   */
  public static void eigen(Complex A[][], Complex lambda[], Complex vec[][], boolean fail[]) {
    //System.out.println("Eigen.eigen(A, lambda, vec, fail)");
    // driver for computing eigenvalues and eigenvectors
    if((A==null)||(lambda==null)||(vec==null)) {
      System.out.println("Error in Eigen.eigen,"+" null or inconsistent array sizes."); //$NON-NLS-1$ //$NON-NLS-2$
      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."); //$NON-NLS-1$ //$NON-NLS-2$
      return;
    }
    fail[0] = false;
    // special cases
    if(n<1) {
      System.out.println("zero size matrix"); //$NON-NLS-1$
      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].mul(A[1][1]);
    c = c.subtract(A[0][1].mul(A[1][0]));
    rad = (b.mul(b)).subtract(c.mul(4.0)); // a==1
    rad = rad.sqrt();
    l1 = (b.add(rad)).div(2.0);
    l2 = (b.subtract(rad)).div(2.0);
    lambda[0] = l1;
    lambda[1] = l2;
    // eigenvectors in columns
    Z[0] = A[0][1].neg();
    Z[1] = A[0][0].subtract(l1);
    t = ComplexMatrix.norm2(Z);
    vec[0][0] = Z[0].div(t);
    vec[1][0] = Z[1].div(t);
    Z[0] = A[1][1].subtract(l2);
    Z[1] = A[1][0].neg();
    t = ComplexMatrix.norm2(Z);
    vec[0][1] = Z[0].div(t);
    vec[1][1] = Z[1].div(t);
  }

  private static double sumabs(Complex Z) {
    return Math.abs(Z.re())+Math.abs(Z.im());
  } // 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.div(x);
            A[ii][m-1] = y;
            for(int j = m; j<n; j++) {
              A[ii][j] = A[ii][j].subtract(y.mul(A[m][j]));
            }
            for(int j = 0; j<n; j++) {
              A[j][m] = A[j][m].add(y.mul(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].re())+Math.abs(A[ien-1][ien-2].re()), Math.abs(A[ien][ien-1].im())+Math.abs(A[ien-1][ien-2].im()));
        } else {
          S = A[ien][ien];
          x = A[ien-1][ien].mul(A[ien][ien-1]);
          if(sumabs(x)>0.0) {
            y = (A[ien-1][ien-1].subtract(S)).div(new Complex(2.0, 0.0));
            z = ((y.mul(y)).add(x)).sqrt();
            if(y.re()*z.re()+y.im()*z.im()<0.0) {
              z = z.neg();
            }
            yy = y.add(z);
            S = S.subtract(x.div(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.div(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.div(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.mul(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].re()>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.mul(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.mul(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].mul(A[jj][ien]));
        }
        y = x.subtract(lambda[i]);
        if(sumabs(y)==0.0) {
          y = accnorm;
        }
        A[i][ien] = z.div(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].mul(A[k][jj]));
        }
        vec[i][jj] = z;
      }
    }
  }                                  // end cxeig2c

}

/*
 * Open Source Physics software is free software; you can redistribute
 * it and/or modify it under the terms of the GNU General Public License (GPL) as
 * published by the Free Software Foundation; either version 2 of the License,
 * or(at your option) any later version.

 * Code that uses any portion of the code in the org.opensourcephysics package
 * or any subpackage (subdirectory) of this package must must also be be released
 * under the GNU GPL license.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
 * or view the license online at http://www.gnu.org/copyleft/gpl.html
 *
 * Copyright (c) 2007  The Open Source Physics project
 *                     http://www.opensourcephysics.org
 */
