content: /*******************************************************************************
   * Copyright (c) 2010 Haifeng Li
   *
   * Licensed under the Apache License, Version 2.0 (the "License");
   * you may not use this file except in compliance with the License.
   * You may obtain a copy of the License at
   *
   *     http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   *******************************************************************************/
  
  package smile.math.matrix;
  
  import smile.math.Complex;
  import smile.math.Math;
  
  import java.util.Arrays;
  
  /**
   * A pure Java implementation of DenseMatrix whose data is stored in a single 1D array of
   * doubles in column major order.
   */
  public class JMatrix implements DenseMatrix {
      private static final long serialVersionUID = 1L;
  
      /**
       * The matrix storage.
       */
      private double[] A;
      /**
       * The number of rows.
       */
      private int nrows;
      /**
       * The number of columns.
       */
      private int ncols;
      /**
       * True if the matrix is symmetric.
       */
      private boolean symmetric = false;
  
      /**
       * Constructor.
       * @param A the array of matrix.
       */
      public JMatrix(double[][] A) {
          this.nrows = A.length;
          this.ncols = A[0].length;
          this.A = new double[nrows*ncols];
  
          for (int i = 0; i < nrows; i++) {
              for (int j = 0; j < ncols; j++) {
                  set(i, j, A[i][j]);
              }
          }
      }
  
      /**
       * Constructor of a column vector/matrix with given array as the internal storage.
       * @param A the array of column vector.
       */
      public JMatrix(double[] A) {
          this.nrows = A.length;
          this.ncols = 1;
          this.A = A;
      }
  
      /**
       * Constructor of all-zero matrix.
       */
      public JMatrix(int rows, int cols) {
          this.nrows = rows;
          this.ncols = cols;
          A = new double[rows*cols];
      }
  
      /**
       * Constructor. Fill the matrix with given value.
       */
      public JMatrix(int rows, int cols, double value) {
          this(rows, cols);
          Arrays.fill(A, value);
      }
  
      /**
       * Constructor.
       * @param value the array of matrix values arranged in column major format
       */
      public JMatrix(int rows, int cols, double[] value) {
          this.nrows = rows;
          this.ncols = cols;
          this.A = value;
      }
  
      @Override
      public boolean isSymmetric() {
          return symmetric;
      }
  
      @Override
      public void setSymmetric(boolean symmetric) {
          this.symmetric = symmetric;
      }
  
      @Override
      public JMatrix copy() {
          JMatrix a = new JMatrix(nrows, ncols, A.clone());
          a.setSymmetric(isSymmetric());
          return a;
      }
  
      @Override
      public double[] data() {
          return A;
      }
  
      @Override
      public JMatrix transpose() {
          JMatrix B = new JMatrix(ncols(), nrows());
          for (int i = 0; i < nrows(); i++) {
              for (int j = 0; j < ncols(); j++) {
                  B.set(j, i, get(i, j));
              }
          }
  
          return B;
      }
  
      @Override
      public String toString() {
          return toString(false);
      }
  
      @Override
      public int nrows() {
          return nrows;
      }
  
      @Override
      public int ncols() {
          return ncols;
      }
  
      @Override
      public int ld() {
          return nrows;
      }
  
      @Override
      public double get(int i, int j) {
          return A[j*nrows + i];
      }
  
      @Override
      public double set(int i, int j, double x) {
          return A[j*nrows + i] = x;
      }
  
      @Override
      public double add(int i, int j, double x) {
          return A[j*nrows + i] += x;
      }
  
      @Override
      public double sub(int i, int j, double x) {
          return A[j*nrows + i] -= x;
      }
  
      @Override
      public double mul(int i, int j, double x) {
          return A[j*nrows + i] *= x;
      }
  
      @Override
      public double div(int i, int j, double x) {
          return A[j*nrows + i] /= x;
      }
  
      @Override
      public JMatrix add(DenseMatrix b) {
          if (b instanceof JMatrix) {
              return add((JMatrix) b);
          } else {
              if (nrows() != b.nrows() || ncols() != b.ncols()) {
                  throw new IllegalArgumentException("Matrix is not of same size.");
              }
  
              int m = nrows();
              int n = ncols();
  
              for (int j = 0; j < n; j++) {
                  for (int i = 0; i < m; i++) {
                      add(i, j, b.get(i, j));
                  }
              }
              return this;
          }
      }
  
      public JMatrix add(JMatrix b) {
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              A[i] += b.A[i];
          }
          return this;
      }
  
      @Override
      public DenseMatrix add(DenseMatrix b, DenseMatrix c) {
          if (b instanceof JMatrix && c instanceof JMatrix) {
              return add((JMatrix) b, (JMatrix) c);
          }
  
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          int m = nrows();
          int n = ncols();
  
          for (int j = 0; j < n; j++) {
              for (int i = 0; i < m; i++) {
                  c.set(i, j, get(i, j) + b.get(i, j));
              }
          }
          return c;
      }
  
      public JMatrix add(JMatrix b, JMatrix c) {
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              c.A[i] = A[i] + b.A[i];
          }
          return c;
      }
  
      @Override
      public JMatrix sub(DenseMatrix b) {
          if (b instanceof JMatrix) {
              return sub((JMatrix) b);
          }
  
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          int m = nrows();
          int n = ncols();
  
          for (int j = 0; j < n; j++) {
              for (int i = 0; i < m; i++) {
                  sub(i, j, b.get(i, j));
              }
          }
          return this;
      }
  
      public JMatrix sub(JMatrix b) {
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              A[i] -= b.A[i];
          }
          return this;
      }
  
  
      @Override
      public DenseMatrix sub(DenseMatrix b, DenseMatrix c) {
          if (b instanceof JMatrix && c instanceof JMatrix) {
              return sub((JMatrix) b, (JMatrix) c);
          }
  
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          int m = nrows();
          int n = ncols();
  
          for (int j = 0; j < n; j++) {
              for (int i = 0; i < m; i++) {
                  c.set(i, j, get(i, j) - b.get(i, j));
              }
          }
          return c;
      }
  
      public JMatrix sub(JMatrix b, JMatrix c) {
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              c.A[i] = A[i] - b.A[i];
          }
          return c;
      }
  
      @Override
      public JMatrix mul(DenseMatrix b) {
          if (b instanceof JMatrix) {
              return mul((JMatrix) b);
          }
  
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          int m = nrows();
          int n = ncols();
  
          for (int j = 0; j < n; j++) {
              for (int i = 0; i < m; i++) {
                  mul(i, j, b.get(i, j));
              }
          }
          return this;
      }
  
      public JMatrix mul(JMatrix b) {
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              A[i] *= b.A[i];
          }
          return this;
      }
  
      @Override
      public DenseMatrix mul(DenseMatrix b, DenseMatrix c) {
          if (b instanceof JMatrix && c instanceof JMatrix) {
              return mul((JMatrix) b, (JMatrix) c);
          }
  
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          int m = nrows();
          int n = ncols();
  
          for (int j = 0; j < n; j++) {
              for (int i = 0; i < m; i++) {
                  c.set(i, j, get(i, j) * b.get(i, j));
              }
          }
          return c;
      }
  
      public JMatrix mul(JMatrix b, JMatrix c) {
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              c.A[i] = A[i] * b.A[i];
          }
          return c;
      }
  
      @Override
      public JMatrix div(DenseMatrix b) {
          if (b instanceof JMatrix) {
              return div((JMatrix) b);
          }
  
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          int m = nrows();
          int n = ncols();
  
          for (int j = 0; j < n; j++) {
              for (int i = 0; i < m; i++) {
                  div(i, j, b.get(i, j));
              }
          }
          return this;
      }
  
      public JMatrix div(JMatrix b) {
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              A[i] /= b.A[i];
          }
          return this;
      }
  
      @Override
      public DenseMatrix div(DenseMatrix b, DenseMatrix c) {
          if (b instanceof JMatrix && c instanceof JMatrix) {
              return div((JMatrix) b, (JMatrix) c);
          }
  
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          int m = nrows();
          int n = ncols();
  
          for (int j = 0; j < n; j++) {
              for (int i = 0; i < m; i++) {
                  c.set(i, j, get(i, j) / b.get(i, j));
              }
          }
          return c;
      }
  
      public JMatrix div(JMatrix b, JMatrix c) {
          if (nrows() != b.nrows() || ncols() != b.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              c.A[i] = A[i] / b.A[i];
          }
          return c;
      }
  
      @Override
      public JMatrix add(double x) {
          for (int i = 0; i < A.length; i++) {
              A[i] += x;
          }
  
          return this;
      }
  
      @Override
      public DenseMatrix add(double x, DenseMatrix c) {
          if (c instanceof JMatrix) {
              return add(x, (JMatrix)c);
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int j = 0; j < ncols; j++) {
              for (int i = 0; i < nrows; i++) {
                  c.set(i, j, get(i, j) + x);
              }
          }
  
          return c;
      }
  
      public JMatrix add(double x, JMatrix c) {
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              c.A[i] = A[i] + x;
          }
  
          return c;
      }
  
      @Override
      public JMatrix sub(double x) {
          for (int i = 0; i < A.length; i++) {
              A[i] -= x;
          }
  
          return this;
      }
  
      @Override
      public DenseMatrix sub(double x, DenseMatrix c) {
          if (c instanceof JMatrix) {
              return sub(x, (JMatrix) c);
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int j = 0; j < ncols; j++) {
              for (int i = 0; i < nrows; i++) {
                  c.set(i, j, get(i, j) - x);
              }
          }
  
          return c;
      }
  
      public JMatrix sub(double x, JMatrix c) {
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              c.A[i] = A[i] - x;
          }
  
          return c;
      }
  
      @Override
      public JMatrix mul(double x) {
          for (int i = 0; i < A.length; i++) {
              A[i] *= x;
          }
  
          return this;
      }
  
      @Override
      public DenseMatrix mul(double x, DenseMatrix c) {
          if (c instanceof JMatrix) {
              return mul(x, (JMatrix)c);
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int j = 0; j < ncols; j++) {
              for (int i = 0; i < nrows; i++) {
                  c.set(i, j, get(i, j) * x);
              }
          }
  
          return c;
      }
  
      public JMatrix mul(double x, JMatrix c) {
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              c.A[i] = A[i] * x;
          }
  
          return c;
      }
  
      @Override
      public JMatrix div(double x) {
          for (int i = 0; i < A.length; i++) {
              A[i] /= x;
          }
  
          return this;
      }
  
      @Override
      public DenseMatrix div(double x, DenseMatrix c) {
          if (c instanceof JMatrix) {
              return div(x, (JMatrix) c);
          }
  
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int j = 0; j < ncols; j++) {
              for (int i = 0; i < nrows; i++) {
                  c.set(i, j, get(i, j) / x);
              }
          }
  
          return c;
      }
  
      public JMatrix div(double x, JMatrix c) {
          if (nrows() != c.nrows() || ncols() != c.ncols()) {
              throw new IllegalArgumentException("Matrix is not of same size.");
          }
  
          for (int i = 0; i < A.length; i++) {
              c.A[i] = A[i] / x;
          }
  
          return c;
      }
  
      @Override
      public JMatrix replaceNaN(double x) {
          for (int i = 0; i < A.length; i++) {
              if (Double.isNaN(A[i])) {
                  A[i] = x;
              }
          }
  
          return this;
      }
  
      @Override
      public double sum() {
          double s = 0.0;
          for (int i = 0; i < A.length; i++) {
              s += A[i];
          }
  
          return s;
      }
  
      @Override
      public JMatrix ata() {
          JMatrix C = new JMatrix(ncols, ncols);
          for (int i = 0; i < ncols; i++) {
              for (int j = 0; j < ncols; j++) {
                  double v = 0.0;
                  for (int k = 0; k < nrows; k++) {
                      v += get(k, i) * get(k, j);
                  }
                  C.set(i, j, v);
              }
          }
          return C;
      }
  
      @Override
      public JMatrix aat() {
          JMatrix C = new JMatrix(nrows, nrows);
          for (int k = 0; k < ncols; k++) {
              for (int i = 0; i < nrows; i++) {
                  for (int j = 0; j < nrows; j++) {
                      C.add(i, j, get(i, k) * get(j, k));
                  }
              }
          }
          return C;
      }
  
      @Override
      public double[] ax(double[] x, double[] y) {
          int n = Math.min(nrows, y.length);
          int p = Math.min(ncols, x.length);
  
          Arrays.fill(y, 0.0);
          for (int k = 0; k < p; k++) {
              for (int i = 0; i < n; i++) {
                  y[i] += get(i, k) * x[k];
              }
          }
  
          return y;
      }
  
      @Override
      public double[] axpy(double[] x, double[] y) {
          int n = Math.min(nrows, y.length);
          int p = Math.min(ncols, x.length);
  
          for (int k = 0; k < p; k++) {
              for (int i = 0; i < n; i++) {
                  y[i] += get(i, k) * x[k];
              }
          }
  
          return y;
      }
  
      @Override
      public double[] axpy(double[] x, double[] y, double b) {
          int n = Math.min(nrows, y.length);
          int p = Math.min(ncols, x.length);
  
          for (int i = 0; i < n; i++) {
              y[i] *= b;
          }
  
          for (int k = 0; k < p; k++) {
              for (int i = 0; i < n; i++) {
                  y[i] += get(i, k) * x[k];
              }
          }
  
          return y;
      }
  
      @Override
      public double[] atx(double[] x, double[] y) {
          int n = Math.min(ncols, y.length);
          int p = Math.min(nrows, x.length);
  
          Arrays.fill(y, 0.0);
          for (int i = 0; i < n; i++) {
              for (int k = 0; k < p; k++) {
                  y[i] += get(k, i) * x[k];
              }
          }
  
          return y;
      }
  
      @Override
      public double[] atxpy(double[] x, double[] y) {
          int n = Math.min(ncols, y.length);
          int p = Math.min(nrows, x.length);
  
          for (int i = 0; i < n; i++) {
              for (int k = 0; k < p; k++) {
                  y[i] += get(k, i) * x[k];
              }
          }
  
          return y;
      }
  
      @Override
      public double[] atxpy(double[] x, double[] y, double b) {
          int n = Math.min(ncols, y.length);
          int p = Math.min(nrows, x.length);
  
          for (int i = 0; i < n; i++) {
              y[i] *= b;
              for (int k = 0; k < p; k++) {
                  y[i] += get(k, i) * x[k];
              }
          }
  
          return y;
      }
  
      @Override
      public JMatrix abmm(DenseMatrix B) {
          if (ncols() != B.nrows()) {
              throw new IllegalArgumentException(String.format("Matrix multiplication A * B: %d x %d vs %d x %d", nrows(), ncols(), B.nrows(), B.ncols()));
          }
  
          JMatrix C = new JMatrix(nrows, B.ncols());
          for (int i = 0; i < nrows; i++) {
              for (int j = 0; j < B.ncols(); j++) {
                  double v = 0.0;
                  for (int k = 0; k < ncols; k++) {
                      v += get(i, k) * B.get(k, j);
                  }
                  C.set(i, j, v);
              }
          }
          return C;
      }
  
      @Override
      public JMatrix abtmm(DenseMatrix B) {
          if (ncols() != B.ncols()) {
              throw new IllegalArgumentException(String.format("Matrix multiplication A * B': %d x %d vs %d x %d", nrows(), ncols(), B.nrows(), B.ncols()));
          }
  
          JMatrix C = new JMatrix(nrows, B.nrows());
          for (int k = 0; k < ncols; k++) {
              for (int i = 0; i < nrows; i++) {
                  for (int j = 0; j < B.nrows(); j++) {
                      C.add(i, j, get(i, k) * B.get(j, k));
                  }
              }
          }
          return C;
      }
  
      @Override
      public JMatrix atbmm(DenseMatrix B) {
          if (nrows() != B.nrows()) {
              throw new IllegalArgumentException(String.format("Matrix multiplication A' * B: %d x %d vs %d x %d", nrows(), ncols(), B.nrows(), B.ncols()));
          }
  
          JMatrix C = new JMatrix(ncols, B.ncols());
          for (int i = 0; i < ncols; i++) {
              for (int j = 0; j < B.ncols(); j++) {
                  double v = 0.0;
                  for (int k = 0; k < nrows; k++) {
                      v += get(k, i) * B.get(k, j);
                  }
                  C.set(i, j, v);
              }
          }
          return C;
      }
  
      /**
       * LU decomposition is computed by a "left-looking", dot-product, Crout/Doolittle algorithm.
       */
      @Override
      public LU lu() {
          int m = nrows();
          int n = ncols();
  
          int[] piv = new int[m];
          for (int i = 0; i < m; i++) {
              piv[i] = i;
          }
  
          int pivsign = 1;
          double[] LUcolj = new double[m];
  
          for (int j = 0; j < n; j++) {
  
              // Make a copy of the j-th column to localize references.
              for (int i = 0; i < m; i++) {
                  LUcolj[i] = get(i, j);
              }
  
              // Apply previous transformations.
              for (int i = 0; i < m; i++) {
                  // Most of the time is spent in the following dot product.
  
                  int kmax = Math.min(i, j);
                  double s = 0.0;
                  for (int k = 0; k < kmax; k++) {
                      s += get(i, k) * LUcolj[k];
                  }
  
                  LUcolj[i] -= s;
                  set(i, j, LUcolj[i]);
              }
  
              // Find pivot and exchange if necessary.
              int p = j;
              for (int i = j + 1; i < m; i++) {
                  if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
                      p = i;
                  }
              }
              if (p != j) {
                  for (int k = 0; k < n; k++) {
                      double t = get(p, k);
                      set(p, k, get(j, k));
                      set(j, k, t);
                  }
                  int k = piv[p];
                  piv[p] = piv[j];
                  piv[j] = k;
                  pivsign = -pivsign;
              }
  
              // Compute multipliers.
              if (j < m & get(j, j) != 0.0) {
                  for (int i = j + 1; i < m; i++) {
                      div(i, j, get(j, j));
                  }
              }
          }
  
          boolean singular = false;
          for (int j = 0; j < n; j++) {
              if (get(j, j) == 0) {
                  singular = true;
                  break;
              }
          }
  
          return new LU(this, piv, pivsign, singular);
      }
  
      /**
       * Cholesky decomposition for symmetric and positive definite matrix.
       * Only the lower triangular part will be used in the decomposition.
       *
       * @throws IllegalArgumentException if the matrix is not positive definite.
       */
      @Override
      public Cholesky cholesky() {
          if (nrows() != ncols()) {
              throw new UnsupportedOperationException("Cholesky decomposition on non-square matrix");
          }
  
          int n = nrows();
  
          // Main loop.
          for (int j = 0; j < n; j++) {
              double d = 0.0;
              for (int k = 0; k < j; k++) {
                  double s = 0.0;
                  for (int i = 0; i < k; i++) {
                      s += get(k, i) * get(j, i);
                  }
                  s = (get(j, k) - s) / get(k, k);
                  set(j, k, s);
                  d = d + s * s;
              }
              d = get(j, j) - d;
  
              if (d < 0.0) {
                  throw new IllegalArgumentException("The matrix is not positive definite.");
              }
  
              set(j, j, Math.sqrt(d));
          }
  
          return new Cholesky(this);
      }
  
      /**
       * QR Decomposition is computed by Householder reflections.
       */
      @Override
      public QR qr() {
          // Initialize.
          int m = nrows();
          int n = ncols();
          double[] rDiagonal = new double[n];
  
          // Main loop.
          for (int k = 0; k < n; k++) {
              // Compute 2-norm of k-th column without under/overflow.
              double nrm = 0.0;
              for (int i = k; i < m; i++) {
                  nrm = Math.hypot(nrm, get(i, k));
              }
  
              if (nrm != 0.0) {
                  // Form k-th Householder vector.
                  if (get(k, k) < 0) {
                      nrm = -nrm;
                  }
                  for (int i = k; i < m; i++) {
                      div(i, k, nrm);
                  }
                  add(k, k, 1.0);
  
                  // Apply transformation to remaining columns.
                  for (int j = k + 1; j < n; j++) {
                      double s = 0.0;
                      for (int i = k; i < m; i++) {
                          s += get(i, k) * get(i, j);
                      }
                      s = -s / get(k, k);
                      for (int i = k; i < m; i++) {
                          add(i, j, s * get(i, k));
                      }
                  }
              }
              rDiagonal[k] = -nrm;
          }
  
          boolean singular = false;
          for (int j = 0; j < rDiagonal.length; j++) {
              if (rDiagonal[j] == 0) {
                  singular = true;
                  break;
              }
          }
  
          return new QR(this, rDiagonal, singular);
      }
  
      @Override
      public SVD svd() {
          int m = nrows();
          int n = ncols();
  
          boolean flag;
          int i, its, j, jj, k, l = 0, nm = 0;
          double anorm, c, f, g, h, s, scale, x, y, z;
          g = scale = anorm = 0.0;
  
          DenseMatrix U = this;
          DenseMatrix V = Matrix.zeros(n, n);
          double[] w = new double[n];
          double[] rv1 = new double[n];
  
          for (i = 0; i < n; i++) {
              l = i + 2;
              rv1[i] = scale * g;
              g = s = scale = 0.0;
  
              if (i < m) {
                  for (k = i; k < m; k++) {
                      scale += Math.abs(U.get(k, i));
                  }
  
                  if (scale != 0.0) {
                      for (k = i; k < m; k++) {
                          U.div(k, i, scale);
                          s += U.get(k, i) * U.get(k, i);
                      }
  
                      f = U.get(i, i);
                      g = -Math.copySign(Math.sqrt(s), f);
                      h = f * g - s;
                      U.set(i, i, f - g);
                      for (j = l - 1; j < n; j++) {
                          for (s = 0.0, k = i; k < m; k++) {
                              s += U.get(k, i) * U.get(k, j);
                          }
                          f = s / h;
                          for (k = i; k < m; k++) {
                              U.add(k, j, f * U.get(k, i));
                          }
                      }
                      for (k = i; k < m; k++) {
                          U.mul(k, i, scale);
                      }
                  }
              }
  
              w[i] = scale * g;
              g = s = scale = 0.0;
  
              if (i + 1 <= m && i + 1 != n) {
                  for (k = l - 1; k < n; k++) {
                      scale += Math.abs(U.get(i, k));
                  }
  
                  if (scale != 0.0) {
                      for (k = l - 1; k < n; k++) {
                          U.div(i, k, scale);
                          s += U.get(i, k) * U.get(i, k);
                      }
  
                      f = U.get(i, l - 1);
                      g = -Math.copySign(Math.sqrt(s), f);
                      h = f * g - s;
                      U.set(i, l - 1, f - g);
  
                      for (k = l - 1; k < n; k++) {
                          rv1[k] = U.get(i, k) / h;
                      }
  
                      for (j = l - 1; j < m; j++) {
                          for (s = 0.0, k = l - 1; k < n; k++) {
                              s += U.get(j, k) * U.get(i, k);
                          }
  
                          for (k = l - 1; k < n; k++) {
                              U.add(j, k, s * rv1[k]);
                          }
                      }
  
                      for (k = l - 1; k < n; k++) {
                          U.mul(i, k, scale);
                      }
                  }
              }
  
              anorm = Math.max(anorm, (Math.abs(w[i]) + Math.abs(rv1[i])));
          }
  
          for (i = n - 1; i >= 0; i--) {
              if (i < n - 1) {
                  if (g != 0.0) {
                      for (j = l; j < n; j++) {
                          V.set(j, i, (U.get(i, j) / U.get(i, l)) / g);
                      }
                      for (j = l; j < n; j++) {
                          for (s = 0.0, k = l; k < n; k++) {
                              s += U.get(i, k) * V.get(k, j);
                          }
                          for (k = l; k < n; k++) {
                              V.add(k, j, s * V.get(k, i));
                          }
                      }
                  }
                  for (j = l; j < n; j++) {
                      V.set(i, j, 0.0);
                      V.set(j, i, 0.0);
                  }
              }
              V.set(i, i, 1.0);
              g = rv1[i];
              l = i;
          }
  
          for (i = Math.min(m, n) - 1; i >= 0; i--) {
              l = i + 1;
              g = w[i];
              for (j = l; j < n; j++) {
                  U.set(i, j, 0.0);
              }
  
              if (g != 0.0) {
                  g = 1.0 / g;
                  for (j = l; j < n; j++) {
                      for (s = 0.0, k = l; k < m; k++) {
                          s += U.get(k, i) * U.get(k, j);
                      }
                      f = (s / U.get(i, i)) * g;
                      for (k = i; k < m; k++) {
                          U.add(k, j, f * U.get(k, i));
                      }
                  }
                  for (j = i; j < m; j++) {
                      U.mul(j, i, g);
                  }
              } else {
                  for (j = i; j < m; j++) {
                      U.set(j, i, 0.0);
                  }
              }
  
              U.add(i, i, 1.0);
          }
  
          for (k = n - 1; k >= 0; k--) {
              for (its = 0; its < 30; its++) {
                  flag = true;
                  for (l = k; l >= 0; l--) {
                      nm = l - 1;
                      if (l == 0 || Math.abs(rv1[l]) <= Math.EPSILON * anorm) {
                          flag = false;
                          break;
                      }
                      if (Math.abs(w[nm]) <= Math.EPSILON * anorm) {
                          break;
                      }
                  }
  
                  if (flag) {
                      c = 0.0;
                      s = 1.0;
                      for (i = l; i < k + 1; i++) {
                          f = s * rv1[i];
                          rv1[i] = c * rv1[i];
                          if (Math.abs(f) <= Math.EPSILON * anorm) {
                              break;
                          }
                          g = w[i];
                          h = Math.hypot(f, g);
                          w[i] = h;
                          h = 1.0 / h;
                          c = g * h;
                          s = -f * h;
                          for (j = 0; j < m; j++) {
                              y = U.get(j, nm);
                              z = U.get(j, i);
                              U.set(j, nm, y * c + z * s);
                              U.set(j,  i, z * c - y * s);
                          }
                      }
                  }
  
                  z = w[k];
                  if (l == k) {
                      if (z < 0.0) {
                          w[k] = -z;
                          for (j = 0; j < n; j++) {
                              V.set(j, k, -V.get(j, k));
                          }
                      }
                      break;
                  }
  
                  if (its == 29) {
                      throw new IllegalStateException("no convergence in 30 iterations");
                  }
  
                  x = w[l];
                  nm = k - 1;
                  y = w[nm];
                  g = rv1[nm];
                  h = rv1[k];
                  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
                  g = Math.hypot(f, 1.0);
                  f = ((x - z) * (x + z) + h * ((y / (f + Math.copySign(g, f))) - h)) / x;
                  c = s = 1.0;
  
                  for (j = l; j <= nm; j++) {
                      i = j + 1;
                      g = rv1[i];
                      y = w[i];
                      h = s * g;
                      g = c * g;
                      z = Math.hypot(f, h);
                      rv1[j] = z;
                      c = f / z;
                      s = h / z;
                      f = x * c + g * s;
                      g = g * c - x * s;
                      h = y * s;
                      y *= c;
  
                      for (jj = 0; jj < n; jj++) {
                          x = V.get(jj, j);
                          z = V.get(jj, i);
                          V.set(jj, j, x * c + z * s);
                          V.set(jj, i, z * c - x * s);
                      }
  
                      z = Math.hypot(f, h);
                      w[j] = z;
                      if (z != 0.0) {
                          z = 1.0 / z;
                          c = f * z;
                          s = h * z;
                      }
  
                      f = c * g + s * y;
                      x = c * y - s * g;
                      for (jj = 0; jj < m; jj++) {
                          y = U.get(jj, j);
                          z = U.get(jj, i);
                          U.set(jj, j, y * c + z * s);
                          U.set(jj, i, z * c - y * s);
                      }
                  }
  
                  rv1[l] = 0.0;
                  rv1[k] = f;
                  w[k] = x;
              }
          }
  
          // order singular values
          int inc = 1;
          double sw;
          double[] su = new double[m], sv = new double[n];
  
          do {
              inc *= 3;
              inc++;
          } while (inc <= n);
  
          do {
              inc /= 3;
              for (i = inc; i < n; i++) {
                  sw = w[i];
                  for (k = 0; k < m; k++) {
                      su[k] = U.get(k, i);
                  }
                  for (k = 0; k < n; k++) {
                      sv[k] = V.get(k, i);
                  }
                  j = i;
                  while (w[j - inc] < sw) {
                      w[j] = w[j - inc];
                      for (k = 0; k < m; k++) {
                          U.set(k, j, U.get(k, j - inc));
                      }
                      for (k = 0; k < n; k++) {
                          V.set(k, j, V.get(k, j - inc));
                      }
                      j -= inc;
                      if (j < inc) {
                          break;
                      }
                  }
                  w[j] = sw;
                  for (k = 0; k < m; k++) {
                      U.set(k, j, su[k]);
                  }
                  for (k = 0; k < n; k++) {
                      V.set(k, j, sv[k]);
                  }
  
              }
          } while (inc > 1);
  
          for (k = 0; k < n; k++) {
              s = 0;
              for (i = 0; i < m; i++) {
                  if (U.get(i, k) < 0.) {
                      s++;
                  }
              }
              for (j = 0; j < n; j++) {
                  if (V.get(j, k) < 0.) {
                      s++;
                  }
              }
              if (s > (m + n) / 2) {
                  for (i = 0; i < m; i++) {
                      U.set(i, k, -U.get(i, k));
                  }
                  for (j = 0; j < n; j++) {
                      V.set(j, k, -V.get(j, k));
                  }
              }
          }
  
          // There are at most Math.min(m, n) nonzero singular values.
          double[] sigma = new double[Math.min(m, n)];
          System.arraycopy(w, 0, sigma, 0, sigma.length);
          return new SVD(U, V, sigma);
      }
  
      @Override
      public double[] eig() {
          if (nrows() != ncols()) {
              throw new IllegalArgumentException(String.format("Matrix is not square: %d x %d", nrows(), ncols()));
          }
  
          int n = nrows();
          double[] d = new double[n];
          double[] e = new double[n];
  
          if (isSymmetric()) {
              // Tridiagonalize.
              tred(this, d, e);
              // Diagonalize.
              tql(d, e, n);
          } else {
              double[] scale = balance(this);
              int[] perm = elmhes(this);
              hqr(this, d, e);
              sort(d, e);
          }
  
          double[] w = new double[2*n];
          System.arraycopy(d, 0, w, 0, n);
          System.arraycopy(e, 0, w, n, n);
          return w;
      }
  
      @Override
      public EVD eigen() {
          if (nrows() != ncols()) {
              throw new IllegalArgumentException(String.format("Matrix is not square: %d x %d", nrows(), ncols()));
          }
  
          int n = nrows();
          double[] d = new double[n];
          double[] e = new double[n];
  
          DenseMatrix V;
          if (isSymmetric()) {
              V = this;
  
              // Tridiagonalize.
              tred2(V, d, e);
              // Diagonalize.
              tql2(V, d, e);
  
          } else {
              double[] scale = balance(this);
              int[] perm = elmhes(this);
  
              V = Matrix.eye(n, n);
  
              eltran(this, V, perm);
              hqr2(this, V, d, e);
              balbak(V, scale);
              sort(d, e, V);
          }
  
          return new EVD(V, d, e);
      }
  
      /**
       * Symmetric Householder reduction to tridiagonal form.
       */
      private static void tred(DenseMatrix V, double[] d, double[] e) {
          int n = V.nrows();
          for (int i = 0; i < n; i++) {
              d[i] = V.get(n - 1, i);
          }
  
          // Householder reduction to tridiagonal form.
          for (int i = n - 1; i > 0; i--) {
  
              // Scale to avoid under/overflow.
              double scale = 0.0;
              double h = 0.0;
              for (int k = 0; k < i; k++) {
                  scale = scale + Math.abs(d[k]);
              }
              if (scale == 0.0) {
                  e[i] = d[i - 1];
                  for (int j = 0; j < i; j++) {
                      d[j] = V.get(i - 1, j);
                      V.set(i, j, 0.0);
                      V.set(j, i, 0.0);
                  }
              } else {
  
                  // Generate Householder vector.
                  for (int k = 0; k < i; k++) {
                      d[k] /= scale;
                      h += d[k] * d[k];
                  }
                  double f = d[i - 1];
                  double g = Math.sqrt(h);
                  if (f > 0) {
                      g = -g;
                  }
                  e[i] = scale * g;
                  h = h - f * g;
                  d[i - 1] = f - g;
                  for (int j = 0; j < i; j++) {
                      e[j] = 0.0;
                  }
  
                  // Apply similarity transformation to remaining columns.
                  for (int j = 0; j < i; j++) {
                      f = d[j];
                      V.set(j, i, f);
                      g = e[j] + V.get(j, j) * f;
                      for (int k = j + 1; k <= i - 1; k++) {
                          g += V.get(k, j) * d[k];
                          e[k] += V.get(k, j) * f;
                      }
                      e[j] = g;
                  }
                  f = 0.0;
                  for (int j = 0; j < i; j++) {
                      e[j] /= h;
                      f += e[j] * d[j];
                  }
                  double hh = f / (h + h);
                  for (int j = 0; j < i; j++) {
                      e[j] -= hh * d[j];
                  }
                  for (int j = 0; j < i; j++) {
                      f = d[j];
                      g = e[j];
                      for (int k = j; k <= i - 1; k++) {
                          V.sub(k, j, (f * e[k] + g * d[k]));
                      }
                      d[j] = V.get(i - 1, j);
                      V.set(i, j, 0.0);
                  }
              }
              d[i] = h;
          }
  
          for (int j = 0; j < n; j++) {
              d[j] = V.get(j, j);
          }
          e[0] = 0.0;
      }
  
      /**
       * Symmetric Householder reduction to tridiagonal form.
       */
      private static void tred2(DenseMatrix V, double[] d, double[] e) {
          int n = V.nrows();
          for (int i = 0; i < n; i++) {
              d[i] = V.get(n - 1, i);
          }
  
          // Householder reduction to tridiagonal form.
          for (int i = n - 1; i > 0; i--) {
              // Scale to avoid under/overflow.
              double scale = 0.0;
              double h = 0.0;
              for (int k = 0; k < i; k++) {
                  scale = scale + Math.abs(d[k]);
              }
              if (scale == 0.0) {
                  e[i] = d[i - 1];
                  for (int j = 0; j < i; j++) {
                      d[j] = V.get(i - 1, j);
                      V.set(i, j, 0.0);
                      V.set(j, i, 0.0);
                  }
              } else {
                  // Generate Householder vector.
                  for (int k = 0; k < i; k++) {
                      d[k] /= scale;
                      h += d[k] * d[k];
                  }
                  double f = d[i - 1];
                  double g = Math.sqrt(h);
                  if (f > 0) {
                      g = -g;
                  }
                  e[i] = scale * g;
                  h = h - f * g;
                  d[i - 1] = f - g;
                  for (int j = 0; j < i; j++) {
                      e[j] = 0.0;
                  }
  
                  // Apply similarity transformation to remaining columns.
                  for (int j = 0; j < i; j++) {
                      f = d[j];
                      V.set(j, i, f);
                      g = e[j] + V.get(j, j) * f;
                      for (int k = j + 1; k <= i - 1; k++) {
                          g += V.get(k, j) * d[k];
                          e[k] += V.get(k, j) * f;
                      }
                      e[j] = g;
                  }
                  f = 0.0;
                  for (int j = 0; j < i; j++) {
                      e[j] /= h;
                      f += e[j] * d[j];
                  }
                  double hh = f / (h + h);
                  for (int j = 0; j < i; j++) {
                      e[j] -= hh * d[j];
                  }
                  for (int j = 0; j < i; j++) {
                      f = d[j];
                      g = e[j];
                      for (int k = j; k <= i - 1; k++) {
                          V.sub(k, j, (f * e[k] + g * d[k]));
                      }
                      d[j] = V.get(i - 1, j);
                      V.set(i, j, 0.0);
                  }
              }
              d[i] = h;
          }
  
          // Accumulate transformations.
          for (int i = 0; i < n - 1; i++) {
              V.set(n - 1, i, V.get(i, i));
              V.set(i, i, 1.0);
              double h = d[i + 1];
              if (h != 0.0) {
                  for (int k = 0; k <= i; k++) {
                      d[k] = V.get(k, i + 1) / h;
                  }
                  for (int j = 0; j <= i; j++) {
                      double g = 0.0;
                      for (int k = 0; k <= i; k++) {
                          g += V.get(k, i + 1) * V.get(k, j);
                      }
                      for (int k = 0; k <= i; k++) {
                          V.sub(k, j,  g * d[k]);
                      }
                  }
              }
              for (int k = 0; k <= i; k++) {
                  V.set(k, i + 1, 0.0);
              }
          }
          for (int j = 0; j < n; j++) {
              d[j] = V.get(n - 1, j);
              V.set(n - 1, j, 0.0);
          }
          V.set(n - 1, n - 1, 1.0);
          e[0] = 0.0;
      }
  
      /**
       * Tridiagonal QL Implicit routine for computing eigenvalues of a symmetric,
       * real, tridiagonal matrix.
       *
       * The routine works extremely well in practice. The number of iterations for the first few
       * eigenvalues might be 4 or 5, say, but meanwhile the off-diagonal elements in the lower right-hand
       * corner have been reduced too. The later eigenvalues are liberated with very little work. The
       * average number of iterations per eigenvalue is typically 1.3 - 1.6. The operation count per
       * iteration is O(n), with a fairly large effective coefficient, say, ~20n. The total operation count
       * for the diagonalization is then ~20n * (1.3 - 1.6)n = ~30n^2. If the eigenvectors are required,
       * there is an additional, much larger, workload of about 3n^3 operations.
       *
       * @param d on input, it contains the diagonal elements of the tridiagonal matrix.
       * On output, it contains the eigenvalues.
       * @param e on input, it contains the subdiagonal elements of the tridiagonal
       * matrix, with e[0] arbitrary. On output, its contents are destroyed.
       * @param n the size of all parameter arrays.
       */
      private static void tql(double[] d, double[] e, int n) {
          for (int i = 1; i < n; i++) {
              e[i - 1] = e[i];
          }
          e[n - 1] = 0.0;
  
          double f = 0.0;
          double tst1 = 0.0;
          for (int l = 0; l < n; l++) {
  
              // Find small subdiagonal element
              tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
              int m = l;
              for (; m < n; m++) {
                  if (Math.abs(e[m]) <= Math.EPSILON * tst1) {
                      break;
                  }
              }
  
              // If m == l, d[l] is an eigenvalue,
              // otherwise, iterate.
              if (m > l) {
                  int iter = 0;
                  do {
                      if (++iter >= 30) {
                          throw new RuntimeException("Too many iterations");
                      }
  
                      // Compute implicit shift
                      double g = d[l];
                      double p = (d[l + 1] - d[l]) / (2.0 * e[l]);
                      double r = Math.hypot(p, 1.0);
                      if (p < 0) {
                          r = -r;
                      }
                      d[l] = e[l] / (p + r);
                      d[l + 1] = e[l] * (p + r);
                      double dl1 = d[l + 1];
                      double h = g - d[l];
                      for (int i = l + 2; i < n; i++) {
                          d[i] -= h;
                      }
                      f = f + h;
  
                      // Implicit QL transformation.
                      p = d[m];
                      double c = 1.0;
                      double c2 = c;
                      double c3 = c;
                      double el1 = e[l + 1];
                      double s = 0.0;
                      double s2 = 0.0;
                      for (int i = m - 1; i >= l; i--) {
                          c3 = c2;
                          c2 = c;
                          s2 = s;
                          g = c * e[i];
                          h = c * p;
                          r = Math.hypot(p, e[i]);
                          e[i + 1] = s * r;
                          s = e[i] / r;
                          c = p / r;
                          p = c * d[i] - s * g;
                          d[i + 1] = h + s * (c * g + s * d[i]);
                      }
                      p = -s * s2 * c3 * el1 * e[l] / dl1;
                      e[l] = s * p;
                      d[l] = c * p;
  
                      // Check for convergence.
                  } while (Math.abs(e[l]) > Math.EPSILON * tst1);
              }
              d[l] = d[l] + f;
              e[l] = 0.0;
          }
  
          // Sort eigenvalues and corresponding vectors.
          for (int i = 0; i < n - 1; i++) {
              int k = i;
              double p = d[i];
              for (int j = i + 1; j < n; j++) {
                  if (d[j] > p) {
                      k = j;
                      p = d[j];
                  }
              }
              if (k != i) {
                  d[k] = d[i];
                  d[i] = p;
              }
          }
      }
  
      /**
       * Tridiagonal QL Implicit routine for computing eigenvalues and eigenvectors of a symmetric,
       * real, tridiagonal matrix.
       *
       * The routine works extremely well in practice. The number of iterations for the first few
       * eigenvalues might be 4 or 5, say, but meanwhile the off-diagonal elements in the lower right-hand
       * corner have been reduced too. The later eigenvalues are liberated with very little work. The
       * average number of iterations per eigenvalue is typically 1.3 - 1.6. The operation count per
       * iteration is O(n), with a fairly large effective coefficient, say, ~20n. The total operation count
       * for the diagonalization is then ~20n * (1.3 - 1.6)n = ~30n^2. If the eigenvectors are required,
       * there is an additional, much larger, workload of about 3n^3 operations.
       *
       * @param V on input, it contains the identity matrix. On output, the kth column
       * of V returns the normalized eigenvector corresponding to d[k].
       * @param d on input, it contains the diagonal elements of the tridiagonal matrix.
       * On output, it contains the eigenvalues.
       * @param e on input, it contains the subdiagonal elements of the tridiagonal
       * matrix, with e[0] arbitrary. On output, its contents are destroyed.
       */
      static void tql2(DenseMatrix V, double[] d, double[] e) {
          int n = V.nrows();
          for (int i = 1; i < n; i++) {
              e[i - 1] = e[i];
          }
          e[n - 1] = 0.0;
  
          double f = 0.0;
          double tst1 = 0.0;
          for (int l = 0; l < n; l++) {
              // Find small subdiagonal element
              tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
              int m = l;
              for (; m < n; m++) {
                  if (Math.abs(e[m]) <= Math.EPSILON * tst1) {
                      break;
                  }
              }
  
              // If m == l, d[l] is an eigenvalue,
              // otherwise, iterate.
              if (m > l) {
                  int iter = 0;
                  do {
                      if (++iter >= 30) {
                          throw new RuntimeException("Too many iterations");
                      }
  
                      // Compute implicit shift
                      double g = d[l];
                      double p = (d[l + 1] - g) / (2.0 * e[l]);
                      double r = Math.hypot(p, 1.0);
                      if (p < 0) {
                          r = -r;
                      }
                      d[l] = e[l] / (p + r);
                      d[l + 1] = e[l] * (p + r);
                      double dl1 = d[l + 1];
                      double h = g - d[l];
                      for (int i = l + 2; i < n; i++) {
                          d[i] -= h;
                      }
                      f = f + h;
  
                      // Implicit QL transformation.
                      p = d[m];
                      double c = 1.0;
                      double c2 = c;
                      double c3 = c;
                      double el1 = e[l + 1];
                      double s = 0.0;
                      double s2 = 0.0;
                      for (int i = m - 1; i >= l; i--) {
                          c3 = c2;
                          c2 = c;
                          s2 = s;
                          g = c * e[i];
                          h = c * p;
                          r = Math.hypot(p, e[i]);
                          e[i + 1] = s * r;
                          s = e[i] / r;
                          c = p / r;
                          p = c * d[i] - s * g;
                          d[i + 1] = h + s * (c * g + s * d[i]);
  
                          // Accumulate transformation.
                          for (int k = 0; k < n; k++) {
                              h = V.get(k, i + 1);
                              V.set(k, i + 1, s * V.get(k, i) + c * h);
                              V.set(k, i,     c * V.get(k, i) - s * h);
                          }
                      }
                      p = -s * s2 * c3 * el1 * e[l] / dl1;
                      e[l] = s * p;
                      d[l] = c * p;
  
                      // Check for convergence.
                  } while (Math.abs(e[l]) > Math.EPSILON * tst1);
              }
              d[l] = d[l] + f;
              e[l] = 0.0;
          }
  
          // Sort eigenvalues and corresponding vectors.
          for (int i = 0; i < n - 1; i++) {
              int k = i;
              double p = d[i];
              for (int j = i + 1; j < n; j++) {
                  if (d[j] > p) {
                      k = j;
                      p = d[j];
                  }
              }
              if (k != i) {
                  d[k] = d[i];
                  d[i] = p;
                  for (int j = 0; j < n; j++) {
                      p = V.get(j, i);
                      V.set(j, i, V.get(j, k));
                      V.set(j, k, p);
                  }
              }
          }
      }
  
      /**
       * Given a square matrix, this routine replaces it by a balanced matrix with
       * identical eigenvalues. A symmetric matrix is already balanced and is
       * unaffected by this procedure.
       */
      private static double[] balance(DenseMatrix A) {
          double sqrdx = Math.RADIX * Math.RADIX;
  
          int n = A.nrows();
  
          double[] scale = new double[n];
          for (int i = 0; i < n; i++) {
              scale[i] = 1.0;
          }
  
          boolean done = false;
          while (!done) {
              done = true;
              for (int i = 0; i < n; i++) {
                  double r = 0.0, c = 0.0;
                  for (int j = 0; j < n; j++) {
                      if (j != i) {
                          c += Math.abs(A.get(j, i));
                          r += Math.abs(A.get(i, j));
                      }
                  }
                  if (c != 0.0 && r != 0.0) {
                      double g = r / Math.RADIX;
                      double f = 1.0;
                      double s = c + r;
                      while (c < g) {
                          f *= Math.RADIX;
                          c *= sqrdx;
                      }
                      g = r * Math.RADIX;
                      while (c > g) {
                          f /= Math.RADIX;
                          c /= sqrdx;
                      }
                      if ((c + r) / f < 0.95 * s) {
                          done = false;
                          g = 1.0 / f;
                          scale[i] *= f;
                          for (int j = 0; j < n; j++) {
                              A.mul(i, j, g);
                          }
                          for (int j = 0; j < n; j++) {
                              A.mul(j, i, f);
                          }
                      }
                  }
              }
          }
  
          return scale;
      }
  
      /**
       * Form the eigenvectors of a real nonsymmetric matrix by back transforming
       * those of the corresponding balanced matrix determined by balance.
       */
      private static void balbak(DenseMatrix V, double[] scale) {
          int n = V.nrows();
          for (int i = 0; i < n; i++) {
              for (int j = 0; j < n; j++) {
                  V.mul(i, j, scale[i]);
              }
          }
      }
  
      /**
       * Reduce a real nonsymmetric matrix to upper Hessenberg form.
       */
      private static int[] elmhes(DenseMatrix A) {
          int n = A.nrows();
          int[] perm = new int[n];
  
          for (int m = 1; m < n - 1; m++) {
              double x = 0.0;
              int i = m;
              for (int j = m; j < n; j++) {
                  if (Math.abs(A.get(j, m - 1)) > Math.abs(x)) {
                      x = A.get(j, m - 1);
                      i = j;
                  }
              }
              perm[m] = i;
              if (i != m) {
                  for (int j = m - 1; j < n; j++) {
                      double swap = A.get(i, j);
                      A.set(i, j, A.get(m, j));
                      A.set(m, j, swap);
                  }
                  for (int j = 0; j < n; j++) {
                      double swap = A.get(j, i);
                      A.set(j, i, A.get(j, m));
                      A.set(j, m, swap);
                  }
              }
              if (x != 0.0) {
                  for (i = m + 1; i < n; i++) {
                      double y = A.get(i, m - 1);
                      if (y != 0.0) {
                          y /= x;
                          A.set(i, m - 1, y);
                          for (int j = m; j < n; j++) {
                              A.sub(i, j, y * A.get(m, j));
                          }
                          for (int j = 0; j < n; j++) {
                              A.add(j, m, y * A.get(j, i));
                          }
                      }
                  }
              }
          }
  
          return perm;
      }
  
      /**
       * Accumulate the stabilized elementary similarity transformations used
       * in the reduction of a real nonsymmetric matrix to upper Hessenberg form by elmhes.
       */
      private static void eltran(DenseMatrix A, DenseMatrix V, int[] perm) {
          int n = A.nrows();
          for (int mp = n - 2; mp > 0; mp--) {
              for (int k = mp + 1; k < n; k++) {
                  V.set(k, mp, A.get(k, mp - 1));
              }
              int i = perm[mp];
              if (i != mp) {
                  for (int j = mp; j < n; j++) {
                      V.set(mp, j, V.get(i, j));
                      V.set(i, j, 0.0);
                  }
                  V.set(i, mp, 1.0);
              }
          }
      }
  
      /**
       * Find all eigenvalues of an upper Hessenberg matrix. On input, A can be
       * exactly as output from elmhes. On output, it is destroyed.
       */
      private static void hqr(DenseMatrix A, double[] d, double[] e) {
          int n = A.nrows();
          int nn, m, l, k, j, its, i, mmin;
          double z, y, x, w, v, u, t, s, r = 0.0, q = 0.0, p = 0.0, anorm = 0.0;
  
          for (i = 0; i < n; i++) {
              for (j = Math.max(i - 1, 0); j < n; j++) {
                  anorm += Math.abs(A.get(i, j));
              }
          }
          nn = n - 1;
          t = 0.0;
          while (nn >= 0) {
              its = 0;
              do {
                  for (l = nn; l > 0; l--) {
                      s = Math.abs(A.get(l - 1, l - 1) + Math.abs(A.get(l, l)));
                      if (s == 0.0) {
                          s = anorm;
                      }
                      if (Math.abs(A.get(l, l - 1)) <= Math.EPSILON * s) {
                          A.set(l, l - 1, 0.0);
                          break;
                      }
                  }
                  x = A.get(nn, nn);
                  if (l == nn) {
                      d[nn--] = x + t;
                  } else {
                      y = A.get(nn - 1, nn - 1);
                      w = A.get(nn, nn - 1) * A.get(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 + Math.copySign(z, p);
                              d[nn - 1] = d[nn] = x + z;
                              if (z != 0.0) {
                                  d[nn] = x - w / z;
                              }
                          } else {
                              d[nn] = x + p;
                              e[nn] = -z;
                              d[nn - 1] = d[nn];
                              e[nn - 1] = -e[nn];
                          }
                          nn -= 2;
                      } else {
                          if (its == 30) {
                              throw new IllegalStateException("Too many iterations in hqr");
                          }
                          if (its == 10 || its == 20) {
                              t += x;
                              for (i = 0; i < nn + 1; i++) {
                                  A.sub(i, i, x);
                              }
                              s = Math.abs(A.get(nn, nn - 1)) + Math.abs(A.get(nn - 1, nn - 2));
                              y = x = 0.75 * s;
                              w = -0.4375 * s * s;
                          }
                          ++its;
                          for (m = nn - 2; m >= l; m--) {
                              z = A.get(m, m);
                              r = x - z;
                              s = y - z;
                              p = (r * s - w) / A.get(m + 1, m) + A.get(m, m + 1);
                              q = A.get(m + 1, m + 1) - z - r - s;
                              r = A.get(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.get(m, m - 1)) * (Math.abs(q) + Math.abs(r));
                              v = Math.abs(p) * (Math.abs(A.get(m - 1, m - 1)) + Math.abs(z) + Math.abs(A.get(m + 1, m + 1)));
                              if (u <= Math.EPSILON * v) {
                                  break;
                              }
                          }
                          for (i = m; i < nn - 1; i++) {
                              A.set(i + 2, i, 0.0);
                              if (i != m) {
                                  A.set(i + 2, i - 1, 0.0);
                              }
                          }
                          for (k = m; k < nn; k++) {
                              if (k != m) {
                                  p = A.get(k, k - 1);
                                  q = A.get(k + 1, k - 1);
                                  r = 0.0;
                                  if (k + 1 != nn) {
                                      r = A.get(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 = Math.copySign(Math.sqrt(p * p + q * q + r * r), p)) != 0.0) {
                                  if (k == m) {
                                      if (l != m) {
                                          A.set(k, k - 1, -A.get(k, k - 1));
                                      }
                                  } else {
                                      A.set(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 + 1; j++) {
                                      p = A.get(k, j) + q * A.get(k + 1, j);
                                      if (k + 1 != nn) {
                                          p += r * A.get(k + 2, j);
                                          A.sub(k + 2, j, p * z);
                                      }
                                      A.sub(k + 1, j, p * y);
                                      A.sub(k, j, p * x);
                                  }
                                  mmin = nn < k + 3 ? nn : k + 3;
                                  for (i = l; i < mmin + 1; i++) {
                                      p = x * A.get(i, k) + y * A.get(i, k + 1);
                                      if (k + 1 != nn) {
                                          p += z * A.get(i, k + 2);
                                          A.sub(i, k + 2, p * r);
                                      }
                                      A.sub(i, k + 1, p * q);
                                      A.sub(i, k, p);
                                  }
                              }
                          }
                      }
                  }
              } while (l + 1 < nn);
          }
      }
  
      /**
       * Finds all eigenvalues of an upper Hessenberg matrix A[0..n-1][0..n-1].
       * On input A can be exactly as output from elmhes and eltran. On output, d and e
       * contain the eigenvalues of A, while V is a matrix whose columns contain
       * the corresponding eigenvectors. The eigenvalues are not sorted, except
       * that complex conjugate pairs appear consecutively with the eigenvalue
       * having the positive imaginary part. For a complex eigenvalue, only the
       * eigenvector corresponding to the eigenvalue with positive imaginary part
       * is stored, with real part in V[0..n-1][i] and imaginary part in V[0..n-1][i+1].
       * The eigenvectors are not normalized.
       */
      private static void hqr2(DenseMatrix A, DenseMatrix V, double[] d, double[] e) {
          int n = A.nrows();
          int nn, m, l, k, j, its, i, mmin, na;
          double z = 0.0, y, x, w, v, u, t, s = 0.0, r = 0.0, q = 0.0, p = 0.0, anorm = 0.0, ra, sa, vr, vi;
  
          for (i = 0; i < n; i++) {
              for (j = Math.max(i - 1, 0); j < n; j++) {
                  anorm += Math.abs(A.get(i, j));
              }
          }
          nn = n - 1;
          t = 0.0;
          while (nn >= 0) {
              its = 0;
              do {
                  for (l = nn; l > 0; l--) {
                      s = Math.abs(A.get(l - 1, l - 1)) + Math.abs(A.get(l, l));
                      if (s == 0.0) {
                          s = anorm;
                      }
                      if (Math.abs(A.get(l, l - 1)) <= Math.EPSILON * s) {
                          A.set(l, l - 1, 0.0);
                          break;
                      }
                  }
                  x = A.get(nn, nn);
                  if (l == nn) {
                      d[nn] = x + t;
                      A.set(nn, nn, x + t);
                      nn--;
                  } else {
                      y = A.get(nn - 1, nn - 1);
                      w = A.get(nn, nn - 1) * A.get(nn - 1, nn);
                      if (l == nn - 1) {
                          p = 0.5 * (y - x);
                          q = p * p + w;
                          z = Math.sqrt(Math.abs(q));
                          x += t;
                          A.set(nn, nn, x );
                          A.set(nn - 1, nn - 1, y + t);
                          if (q >= 0.0) {
                              z = p + Math.copySign(z, p);
                              d[nn - 1] = d[nn] = x + z;
                              if (z != 0.0) {
                                  d[nn] = x - w / z;
                              }
                              x = A.get(nn, nn - 1);
                              s = Math.abs(x) + Math.abs(z);
                              p = x / s;
                              q = z / s;
                              r = Math.sqrt(p * p + q * q);
                              p /= r;
                              q /= r;
                              for (j = nn - 1; j < n; j++) {
                                  z = A.get(nn - 1, j);
                                  A.set(nn - 1, j, q * z + p * A.get(nn, j));
                                  A.set(nn, j, q * A.get(nn, j) - p * z);
                              }
                              for (i = 0; i <= nn; i++) {
                                  z = A.get(i, nn - 1);
                                  A.set(i, nn - 1, q * z + p * A.get(i, nn));
                                  A.set(i, nn, q * A.get(i, nn) - p * z);
                              }
                              for (i = 0; i < n; i++) {
                                  z = V.get(i, nn - 1);
                                  V.set(i, nn - 1, q * z + p * V.get(i, nn));
                                  V.set(i, nn, q * V.get(i, nn) - p * z);
                              }
                          } else {
                              d[nn] = x + p;
                              e[nn] = -z;
                              d[nn - 1] = d[nn];
                              e[nn - 1] = -e[nn];
                          }
                          nn -= 2;
                      } else {
                          if (its == 30) {
                              throw new IllegalArgumentException("Too many iterations in hqr");
                          }
                          if (its == 10 || its == 20) {
                              t += x;
                              for (i = 0; i < nn + 1; i++) {
                                  A.sub(i, i, x);
                              }
                              s = Math.abs(A.get(nn, nn - 1)) + Math.abs(A.get(nn - 1, nn - 2));
                              y = x = 0.75 * s;
                              w = -0.4375 * s * s;
                          }
                          ++its;
                          for (m = nn - 2; m >= l; m--) {
                              z = A.get(m, m);
                              r = x - z;
                              s = y - z;
                              p = (r * s - w) / A.get(m + 1, m) + A.get(m, m + 1);
                              q = A.get(m + 1, m + 1) - z - r - s;
                              r = A.get(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.get(m, m - 1)) * (Math.abs(q) + Math.abs(r));
                              v = Math.abs(p) * (Math.abs(A.get(m - 1, m - 1)) + Math.abs(z) + Math.abs(A.get(m + 1, m + 1)));
                              if (u <= Math.EPSILON * v) {
                                  break;
                              }
                          }
                          for (i = m; i < nn - 1; i++) {
                              A.set(i + 2, i , 0.0);
                              if (i != m) {
                                  A.set(i + 2, i - 1, 0.0);
                              }
                          }
                          for (k = m; k < nn; k++) {
                              if (k != m) {
                                  p = A.get(k, k - 1);
                                  q = A.get(k + 1, k - 1);
                                  r = 0.0;
                                  if (k + 1 != nn) {
                                      r = A.get(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 = Math.copySign(Math.sqrt(p * p + q * q + r * r), p)) != 0.0) {
                                  if (k == m) {
                                      if (l != m) {
                                          A.set(k, k - 1, -A.get(k, k - 1));
                                      }
                                  } else {
                                      A.set(k, k - 1, -s * x);
                                  }
                                  p += s;
                                  x = p / s;
                                  y = q / s;
                                  z = r / s;
                                  q /= p;
                                  r /= p;
                                  for (j = k; j < n; j++) {
                                      p = A.get(k, j) + q * A.get(k + 1, j);
                                      if (k + 1 != nn) {
                                          p += r * A.get(k + 2, j);
                                          A.sub(k + 2, j, p * z);
                                      }
                                      A.sub(k + 1, j, p * y);
                                      A.sub(k, j, p * x);
                                  }
                                  mmin = nn < k + 3 ? nn : k + 3;
                                  for (i = 0; i < mmin + 1; i++) {
                                      p = x * A.get(i, k) + y * A.get(i, k + 1);
                                      if (k + 1 != nn) {
                                          p += z * A.get(i, k + 2);
                                          A.sub(i, k + 2, p * r);
                                      }
                                      A.sub(i, k + 1, p * q);
                                      A.sub(i, k, p);
                                  }
                                  for (i = 0; i < n; i++) {
                                      p = x * V.get(i, k) + y * V.get(i, k + 1);
                                      if (k + 1 != nn) {
                                          p += z * V.get(i, k + 2);
                                          V.sub(i, k + 2, p * r);
                                      }
                                      V.sub(i, k + 1, p * q);
                                      V.sub(i, k, p);
                                  }
                              }
                          }
                      }
                  }
              } while (l + 1 < nn);
          }
  
          if (anorm != 0.0) {
              for (nn = n - 1; nn >= 0; nn--) {
                  p = d[nn];
                  q = e[nn];
                  na = nn - 1;
                  if (q == 0.0) {
                      m = nn;
                      A.set(nn, nn, 1.0);
                      for (i = nn - 1; i >= 0; i--) {
                          w = A.get(i, i) - p;
                          r = 0.0;
                          for (j = m; j <= nn; j++) {
                              r += A.get(i, j) * A.get(j, nn);
                          }
                          if (e[i] < 0.0) {
                              z = w;
                              s = r;
                          } else {
                              m = i;
  
                              if (e[i] == 0.0) {
                                  t = w;
                                  if (t == 0.0) {
                                      t = Math.EPSILON * anorm;
                                  }
                                  A.set(i, nn, -r / t);
                              } else {
                                  x = A.get(i, i + 1);
                                  y = A.get(i + 1, i);
                                  q = Math.sqr(d[i] - p) + Math.sqr(e[i]);
                                  t = (x * s - z * r) / q;
                                  A.set(i, nn, t);
                                  if (Math.abs(x) > Math.abs(z)) {
                                      A.set(i + 1, nn, (-r - w * t) / x);
                                  } else {
                                      A.set(i + 1, nn, (-s - y * t) / z);
                                  }
                              }
                              t = Math.abs(A.get(i, nn));
                              if (Math.EPSILON * t * t > 1) {
                                  for (j = i; j <= nn; j++) {
                                      A.div(j, nn, t);
                                  }
                              }
                          }
                      }
                  } else if (q < 0.0) {
                      m = na;
                      if (Math.abs(A.get(nn, na)) > Math.abs(A.get(na, nn))) {
                          A.set(na, na, q / A.get(nn, na));
                          A.set(na, nn, -(A.get(nn, nn) - p) / A.get(nn, na));
                      } else {
                          Complex temp = cdiv(0.0, -A.get(na, nn), A.get(na, na) - p, q);
                          A.set(na, na, temp.re());
                          A.set(na, nn, temp.im());
                      }
                      A.set(nn, na, 0.0);
                      A.set(nn, nn, 1.0);
                      for (i = nn - 2; i >= 0; i--) {
                          w = A.get(i, i) - p;
                          ra = sa = 0.0;
                          for (j = m; j <= nn; j++) {
                              ra += A.get(i, j) * A.get(j, na);
                              sa += A.get(i, j) * A.get(j, nn);
                          }
                          if (e[i] < 0.0) {
                              z = w;
                              r = ra;
                              s = sa;
                          } else {
                              m = i;
                              if (e[i] == 0.0) {
                                  Complex temp = cdiv(-ra, -sa, w, q);
                                  A.set(i, na, temp.re());
                                  A.set(i, nn, temp.im());
                              } else {
                                  x = A.get(i, i + 1);
                                  y = A.get(i + 1, i);
                                  vr = Math.sqr(d[i] - p) + Math.sqr(e[i]) - q * q;
                                  vi = 2.0 * q * (d[i] - p);
                                  if (vr == 0.0 && vi == 0.0) {
                                      vr = Math.EPSILON * anorm * (Math.abs(w) + Math.abs(q) + Math.abs(x) + Math.abs(y) + Math.abs(z));
                                  }
                                  Complex temp = cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
                                  A.set(i, na, temp.re());
                                  A.set(i, nn, temp.im());
                                  if (Math.abs(x) > Math.abs(z) + Math.abs(q)) {
                                      A.set(i + 1, na, (-ra - w * A.get(i, na) + q * A.get(i, nn)) / x);
                                      A.set(i + 1, nn, (-sa - w * A.get(i, nn) - q * A.get(i, na)) / x);
                                  } else {
                                      temp = cdiv(-r - y * A.get(i, na), -s - y * A.get(i, nn), z, q);
                                      A.set(i + 1, na, temp.re());
                                      A.set(i + 1, nn, temp.im());
                                  }
                              }
                          }
                          t = Math.max(Math.abs(A.get(i, na)), Math.abs(A.get(i, nn)));
                          if (Math.EPSILON * t * t > 1) {
                              for (j = i; j <= nn; j++) {
                                  A.div(j, na, t);
                                  A.div(j, nn, t);
                              }
                          }
                      }
                  }
              }
  
              for (j = n - 1; j >= 0; j--) {
                  for (i = 0; i < n; i++) {
                      z = 0.0;
                      for (k = 0; k <= j; k++) {
                          z += V.get(i, k) * A.get(k, j);
                      }
                      V.set(i, j, z);
                  }
              }
          }
      }
  
      /**
       *  Complex scalar division.
       */
      private static Complex cdiv(double xr, double xi, double yr, double yi) {
          double cdivr, cdivi;
          double r, d;
          if (Math.abs(yr) > Math.abs(yi)) {
              r = yi / yr;
              d = yr + r * yi;
              cdivr = (xr + r * xi) / d;
              cdivi = (xi - r * xr) / d;
          } else {
              r = yr / yi;
              d = yi + r * yr;
              cdivr = (r * xr + xi) / d;
              cdivi = (r * xi - xr) / d;
          }
  
          return new Complex(cdivr, cdivi);
      }
  
      /**
       * Sort eigenvalues.
       */
      protected static void sort(double[] d, double[] e) {
          int i = 0;
          int n = d.length;
          for (int j = 1; j < n; j++) {
              double real = d[j];
              double img = e[j];
              for (i = j - 1; i >= 0; i--) {
                  if (d[i] >= d[j]) {
                      break;
                  }
                  d[i + 1] = d[i];
                  e[i + 1] = e[i];
              }
              d[i + 1] = real;
              e[i + 1] = img;
          }
      }
  
      /**
       * Sort eigenvalues and eigenvectors.
       */
      protected static void sort(double[] d, double[] e, DenseMatrix V) {
          int i = 0;
          int n = d.length;
          double[] temp = new double[n];
          for (int j = 1; j < n; j++) {
              double real = d[j];
              double img = e[j];
              for (int k = 0; k < n; k++) {
                  temp[k] = V.get(k, j);
              }
              for (i = j - 1; i >= 0; i--) {
                  if (d[i] >= d[j]) {
                      break;
                  }
                  d[i + 1] = d[i];
                  e[i + 1] = e[i];
                  for (int k = 0; k < n; k++) {
                      V.set(k, i + 1, V.get(k, i));
                  }
              }
              d[i + 1] = real;
              e[i + 1] = img;
              for (int k = 0; k < n; k++) {
                  V.set(k, i + 1, temp[k]);
              }
          }
      }
  }
  