content: package Jampack;
  
  /**
   * Eig implements the eigenvalue-vector decomposition of
   * of a square matrix.  Specifically given a diagonalizable
   * matrix A, there is a matrix nonsingular matrix X such that
   * <pre>
   *      D = X<sup>-1</sup> AX
   * </pre>
   * is diagonal.  The columns of X are eigenvectors of A corresponding
   * to the diagonal elements of D.  Eig implements X as a Zmat and
   * D as a Zdiagmat.
   * <p/>
   * Warning: if A is defective rounding error will allow Eig to
   * compute a set of eigevectors.  However, the matrix X will
   * be ill conditioned.
   *
   * @author G. W. Stewart
   * @version Pre-alpha
   */
  
  public class Eig {
  
      /**
       * The matrix of eigevectors
       */
      public Zmat X;
  
      /**
       * The diagonal matrix of eigenvalues
       */
      public Zdiagmat D;
  
      /**
       * Creates an eigenvalue-vector decomposition of a square matrix A.
       *
       * @param A The matrix whose decomposition is to be
       *          computed
       * @throws JampackException Thrown if A is not square. <br>
       *                          Passed from below.
       */
  
      public Eig(Zmat A)
              throws JampackException {
  
          int i, j, k;
          double norm, scale;
          Z z, d;
  
          A.getProperties();
  
          if (A.nr != A.nc) {
              throw new JampackException
                      ("Matrix not square.");
          }
  
          int n = A.nr;
  
       /* Compute the Schur decomposition of $A$ and set up T and D. */
  
          Schur S = new Schur(A);
  
  
          Zmat T = S.T;
  
          D = new Zdiagmat(T);
  
          norm = Norm.fro(A);
  
          X = new Zmat(n, n);
  
        /* Compute the eigevectors of T */
  
          for (k = n - 1; k >= 0; k--) {
  
              d = T.get0(k, k);
  
              X.re[k][k] = 1.0;
              X.im[k][k] = 0.0;
  
              for (i = k - 1; i >= 0; i--) {
  
                  X.re[i][k] = -T.re[i][k];
                  X.im[i][k] = -T.im[i][k];
  
                  for (j = i + 1; j < k; j++) {
  
                      X.re[i][k] = X.re[i][k] - T.re[i][j] * X.re[j][k]
                              + T.im[i][j] * X.im[j][k];
                      X.im[i][k] = X.im[i][k] - T.re[i][j] * X.im[j][k]
                              - T.im[i][j] * X.re[j][k];
                  }
  
                  z = T.get0(i, i);
                  z.Minus(z, d);
                  if (z.re == 0.0 && z.im == 0.0) { // perturb zero diagonal
                      z.re = 1.0e-16 * norm;      // to avoid division by zero
                  }
                  z.Div(X.get0(i, k), z);
                  X.put0(i, k, z);
              }
  
          /* Scale the vector so its norm is one. */
  
              scale = 1.0 / Norm.fro(X, X.bx, X.rx, X.bx + k, X.bx + k);
              for (i = 0; i < X.nr; i++) {
                  X.re[i][k] = scale * X.re[i][k];
                  X.im[i][k] = scale * X.im[i][k];
              }
          }
          X = Times.o(S.U, X);
      }
  }
  
  
  