smile.math.matrix.JMatrix Maven / Gradle / Ivy
/*******************************************************************************
* 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]);
}
}
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy