Commit 896c849f authored by Alexis  CRISCUOLO's avatar Alexis CRISCUOLO
Browse files

1.992

parent f2f83fbf
......@@ -34,12 +34,14 @@
package bmge;
import java.util.Arrays;
import java.util.stream.IntStream;
/*
NOTE:
This class is a complete rewriting of the class EigenvalueDecomposition from the package JAMA.
+ JAMA: https://math.nist.gov/javanumerics/jama/
+ EigenvalueDecomposition: https://math.nist.gov/javanumerics/jama/doc/Jama/EigenvalueDecomposition.html
+ JAMA: https://math.nist.gov/javanumerics/jama/
+ EigenvalueDecomposition: https://math.nist.gov/javanumerics/jama/doc/Jama/EigenvalueDecomposition.html
+ EigenvalueDecomposition src: https://github.com/fiji/Jama/blob/master/src/main/java/Jama/EigenvalueDecomposition.java
*/
public class Eigen {
......@@ -50,6 +52,9 @@ public class Eigen {
//## data
private int n; //## NOTE: size
private int n_1; //## NOTE: size-1
private int[] id0n; //## NOTE: int array from 0 to n-1
private int[] id0n_1; //## NOTE: int array from 0 to n-2
private int[] id1n; //## NOTE: int array from 1 to n-1
private boolean issymmetric; //## NOTE: true if symmetric
private double[] d, e; //## NOTE: eigenvalues
private double[][] V; //## NOTE: eigenvectors
......@@ -60,12 +65,16 @@ public class Eigen {
private double c, c2, c3, f, g, h, p, q, r, s, s2, t, w, x, y, z;
private double dl1, el1, exshift, hh, norm, scale, tst1, ra, sa, vr, vi;
private boolean notlast;
private int[] id0i, id0i1;;
private double[][] H;
private double[] ort;
public Eigen (final double[][] mat) {
n_1 = n = mat.length;
--n_1;
id0n = IntStream.range(0, n).toArray();
id0n_1 = IntStream.range(0, n_1).toArray();
id1n = IntStream.range(1, n).toArray();
V = new double[n][n];
inV = new double[n][n];
d = new double[n];
......@@ -73,14 +82,14 @@ public class Eigen {
issymmetric = true;
i = n; while ( issymmetric && (j=--i) >= 0 ) while ( issymmetric && --j >= 0 ) issymmetric = ( mat[i][j] == mat[j][i] );
if ( issymmetric ) {
i = n; while ( --i >= 0 ) System.arraycopy(mat[i], 0, V[i], 0, n);
for (int i: id0n) System.arraycopy(mat[i], 0, V[i], 0, n);
tred2(); //## tridiagonalizing
tql2(); //## diagonalizing
}
else {
H = new double[n][n];
ort = new double[n];
i = n; while ( --i >= 0 ) System.arraycopy(mat[i], 0, H[i], 0, n);
for (int i: id0n) System.arraycopy(mat[i], 0, H[i], 0, n);
orthes(); //## reducing to Hessenberg form.
hqr2(); //## reducing Hessenberg to real Schur form.
}
......@@ -104,38 +113,39 @@ public class Eigen {
//## Symmetric Householder reduction to tridiagonal form.
private void tred2 () {
j = n; while ( --j >= 0 ) d[j] = V[n_1][j];
for (int j: id0n) d[j] = V[n_1][j];
i = n;
while ( --i > 0 ) {
// Scale to avoid under/overflow.
i_1 = i; --i_1; scale = h = 0;
k = i; while ( --k >= 0 ) scale += Math.abs(d[k]);
i_1 = i-1;
id0i = IntStream.range(0, i_1).toArray();
scale = h = 0;
for (int k: id0i) scale += Math.abs(d[k]);
if ( scale == 0.0 ) {
e[i] = d[i_1];
j = i; while ( --j >= 0 ) { d[j] = V[i_1][j]; V[i][j] = 0.0; V[j][i] = 0.0; }
for (int j: id0i) { d[j] = V[i_1][j]; V[i][j] = 0.0; V[j][i] = 0.0; }
}
else {
// Generate Householder vector.
k = i; while ( --k >= 0 ) { d[k] /= scale; h += square(d[k]); }
for (int k: id0i) { d[k] /= scale; h += square(d[k]); }
f = d[i_1];
g = ( f > 0 ) ? -Math.sqrt(h) : Math.sqrt(h);
e[i] = scale * g;
h -= f * g;
d[i_1] = f - g;
j = i; while ( --j >= 0 ) e[j] = 0.0;
for (int j: id0i) e[j] = 0.0;
// Apply similarity transformation to remaining columns.
j = -1;
while ( ++j < i ) {
for (int j: id0i) {
V[j][i] = f = d[j];
g = e[j] + V[j][j] * f;
k = j; while ( ++k <= i_1 ) { g += V[k][j] * d[k]; e[k] += V[k][j] * f; }
for (int k: IntStream.range(j+1, i_1).toArray()) { g += V[k][j] * d[k]; e[k] += V[k][j] * f; }
e[j] = g;
}
f = 0; j = i; while ( --j >= 0 ) { e[j] /= h; f += e[j] * d[j]; }
f = 0;
for (int j: id0i) { e[j] /= h; f += e[j] * d[j]; }
hh = f / (2*h);
j = i; while ( --j >= 0 ) e[j] -= hh * d[j];
j = -1;
while ( ++j < i ) {
for (int j: id0i) e[j] -= hh * d[j];
for (int j: id0i) {
f = d[j]; g = e[j];
k = i; while ( --k >= j ) V[k][j] -= (f * e[k] + g * d[k]);
d[j] = V[i_1][j]; V[i][j] = 0;
......@@ -144,34 +154,32 @@ public class Eigen {
d[i] = h;
}
// Accumulate transformations
i = n_1;
while ( --i >= 0 ) {
for (int i: id0n_1) {
V[n_1][i] = V[i][i];
V[i][i] = 1;
i1 = i+1;
id0i1 = IntStream.range(0, i1).toArray();
if ( (h=d[i1]) != 0.0 ) {
k = i1; while ( --k >= 0 ) d[k] = V[k][i1] / h;
j = i1;
while ( --j >= 0 ) {
for (int k: id0i1) d[k] = V[k][i1] / h;
for (int j: id0i1) {
g = 0;
k = i1; while ( --k >= 0 ) g += V[k][i1] * V[k][j];
k = i1; while ( --k >= 0 ) V[k][j] -= g * d[k];
for (int k: id0i1) g += V[k][i1] * V[k][j];
for (int k: id0i1) V[k][j] -= g * d[k];
}
}
k = i1; while ( --k >= 0 ) V[k][i1] = 0;
for (int k: id0i1) V[k][i1] = 0;
}
j = n; while ( --j >= 0 ) { d[j] = V[n_1][j]; V[n_1][j] = 0; }
for (int j: id0n) { d[j] = V[n_1][j]; V[n_1][j] = 0; }
V[n_1][n_1] = 1;
e[0] = 0.0;
}
//## Symmetric tridiagonal QL algorithm.
private void tql2 () {
i = n; while ( --i > 0 ) e[i-1] = e[i];
for (int i: id1n) e[i-1] = e[i];
e[n_1] = 0.0;
f = tst1 = 0;
l = n;
while ( --l >= 0 ) {
for (int l: id0n) {
// Find small subdiagonal element
tst1 = Math.max(tst1, Math.abs(d[l])+Math.abs(e[l]));
m = l; while ( m < n && Math.abs(e[m]) > EPS*tst1 ) ++m;
......@@ -211,7 +219,7 @@ public class Eigen {
d[i+1] = h + s * (c * g + s * d[i]);
// Accumulate transformation
i1 = i+1;
k = n; while ( --k >= 0 ) { h = V[k][i1]; V[k][i1] = s * V[k][i] + c * h; V[k][i] = c * V[k][i] - s * h; }
for (int k: id0n) { h = V[k][i1]; V[k][i1] = s * V[k][i] + c * h; V[k][i] = c * V[k][i] - s * h; }
}
p = -s * s2 * c3 * el1 * e[l] / dl1;
e[l] = s * p;
......@@ -222,14 +230,13 @@ public class Eigen {
e[l] = 0;
}
// Sort eigenvalues and corresponding vectors
i = n;
while ( --i >= 0 ) {
for (int i: id0n_1) {
p = d[j=k=i];
while ( ++j < n ) if ( d[j] < p ) { k = j; p = d[j]; }
if ( k != i ) {
d[k] = d[i];
d[i] = p;
j = n; while ( --j >= 0 ) { p = V[j][i]; V[j][i] = V[j][k]; V[j][k] = p; }
for (int j: id0n) { p = V[j][i]; V[j][i] = V[j][k]; V[j][k] = p; }
}
}
}
......@@ -249,30 +256,30 @@ public class Eigen {
h -= ort[m] * g;
ort[m] -= g;
// Apply Householder similarity transformation: H = (I-u*u'/h)*H*(I-u*u')/h)
j = m; --j;
j = m-1;
while ( ++j < n ) {
f = 0.0; i = high; ++i; while ( --i >= m ) f += ort[i]*H[i][j];
f /= h; i = m; --i; while ( ++i <= high ) H[i][j] -= f*ort[i];
f = 0.0; i = high1; while ( --i >= m ) f += ort[i]*H[i][j];
f /= h; i = m-1; while ( ++i <= high ) H[i][j] -= f*ort[i];
}
i = high1;
while ( --i >= 0 ) {
f = 0.0; j = high1; while ( --j >= m ) f += ort[j]*H[i][j];
f /= h; j = high1; while ( --j >= m ) H[i][j] -= f*ort[j];
i = -1;
while ( ++i <= high ) {
f = 0.0; j = high1; while ( --j >= m ) f += ort[j]*H[i][j];
f /= h; j = m-1; while ( ++j <= high ) H[i][j] -= f*ort[j];
}
ort[m] = scale * ort[m];
H[m][m-1] = scale*g;
}
}
// Accumulate transformations (Algol's ortran).
i = n; while ( --i >= 0 ) V[i][i] = 1.0;
for (int i: id0n) { Arrays.fill(V[i], 0); V[i][i] = 1.0; }
m = high;
while ( --m >= low1 )
if ( H[m][m-1] != 0.0 ) {
i = m; while ( ++i <= high ) ort[i] = H[i][m-1];
j = m; --j;
j = m-1;
while ( ++j <= high ) {
g = 0; i = high1; while ( --i >= m ) g += ort[i] * V[i][j];
g = (g / ort[m]) / H[m][m-1]; i = high1; while ( --i >= m ) V[i][j] += g * ort[i];
g = 0; i = m-1; while ( ++i <= high ) g += ort[i] * V[i][j];
g = (g / ort[m]) / H[m][m-1]; i = m-1; while ( ++i <= high ) V[i][j] += g * ort[i];
}
}
}
......@@ -304,8 +311,7 @@ public class Eigen {
exshift = p = q = r = s = z = 0;
// Store roots isolated by balanc and compute matrix norm
norm = 0;
i = nn;
while ( --i >= 0 ) {
for (int i: id0n) {
if ( i < low || i > high ) { d[i] = H[i][i]; e[i] = 0.0; }
j = ( i-1 > 0 ) ? i-2 : -1; while ( ++j < nn ) norm += Math.abs(H[i][j]);
}
......@@ -330,14 +336,14 @@ public class Eigen {
z = ( p >= 0 ) ? p + z : p - z; d[n-1] = x + z; d[n] = ( z != 0.0 ) ? x - w/z : d[n-1];
e[n-1] = e[n] = 0; x = H[n][n-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;
j = n; --j; --j; while ( ++j < nn ) { z = H[n-1][j]; H[n-1][j] = q*z + p*H[n][j]; H[n][j] = q*H[n][j] - p*z; } // Row modification
i = n; ++i; while ( --i >= 0 ) { z = H[i][n-1]; H[i][n-1] = q*z + p*H[i][n]; H[i][n] = q*H[i][n] - p*z; } // Column modification
i = low; --i; while ( ++i <= high ) { z = V[i][n-1]; V[i][n-1] = q*z + p*V[i][n]; V[i][n] = q*V[i][n] - p*z; } // Accumulate transformations
j = n-2; while ( ++j < nn ) { z = H[n-1][j]; H[n-1][j] = q*z + p*H[n][j]; H[n][j] = q*H[n][j] - p*z; } // Row modification
i = -1; while ( ++i <= n ) { z = H[i][n-1]; H[i][n-1] = q*z + p*H[i][n]; H[i][n] = q*H[i][n] - p*z; } // Column modification
i = low-1; while ( ++i <= high ) { z = V[i][n-1]; V[i][n-1] = q*z + p*V[i][n]; V[i][n] = q*V[i][n] - p*z; } // Accumulate transformations
}
else { // Complex pair
d[n-1] = x+p; d[n] = x+p; e[n-1] = z; e[n] = -z;
}
--n; --n; //############### WARNING!!
n -= 2; //############### WARNING!!
iter = 0;
}
else { // No convergence yet
......@@ -345,13 +351,13 @@ public class Eigen {
x = H[n][n]; y = ( l < n ) ? H[n-1][n-1] : 0; w = ( l < n ) ? H[n][n-1]*H[n-1][n] : 0;
if ( iter == 10 ) { // Wilkinson's original ad hoc shift
exshift += x;
i = low; --i; while ( ++i <= n ) H[i][i] -= x;
i = low-1; while ( ++i <= n ) H[i][i] -= x;
s = Math.abs(H[n][n-1]) + Math.abs(H[n-1][n-2]);
x = y = 0.75*s; w = -0.4375*s*s;
}
if ( iter == 30 && (s=(square((y-x)/2.0) + w)) > 0 ) { // MATLAB's new ad hoc shift
s = ( y < x ) ? -Math.sqrt(s) : Math.sqrt(s); s = x - w/((y-x)/2.0+s);
i = low; --i; while ( ++i <= n ) H[i][i] -= s;
i = low-1; while ( ++i <= n ) H[i][i] -= s;
exshift += s; x = y = w = 0.964;
}
++iter;
......@@ -364,9 +370,9 @@ public class Eigen {
if ( Math.abs(H[m][m-1])*(Math.abs(q)+Math.abs(r)) < EPS*Math.abs(p)*(Math.abs(H[m-1][m-1])+Math.abs(z)+Math.abs(H[m+1][m+1])) ) break;
--m;
}
i = m; ++i; while ( ++i <= n ) { H[i][i-2] = 0.0; if ( i > m+2 ) H[i][i-3] = 0.0; }
i = m+1; while ( ++i <= n ) { H[i][i-2] = 0.0; if ( i > m+2 ) H[i][i-3] = 0.0; }
// Double QR step involving rows l:n and columns m:n
k = m; --k;
k = m-1;
while ( ++k < n ) {
notlast = ( k != n-1 );
if ( k != m ) {
......@@ -379,21 +385,21 @@ public class Eigen {
if ( k != m ) H[k][k-1] = -s*x; else if (l != m) H[k][k-1] = -H[k][k-1];
p += s; x = p/s; y = q/s; z = r/s; q /= p; r /= p;
// Row modification
j = nn;
while ( --j >= k ) {
j = k-1;
while ( ++j < nn ) {
p = H[k][j] + q*H[k+1][j];
if ( notlast ) { p += r * H[k+2][j]; H[k+2][j] -= p*z; }
H[k][j] -= p*x; H[k+1][j] -= p*y;
}
// Column modification
i = ( n < k+3 ) ? n+1 : k+4;
while ( --i >= 0 ) {
i = -1;
while ( ++i <= ((n<k+3)?n:k+3) ) {
p = x*H[i][k] + y*H[i][k+1];
if ( notlast ) { p += z*H[i][k+2]; H[i][k+2] -= p*r; }
H[i][k] -= p; H[i][k+1] -= p*q;
}
// Accumulate transformations
i = low; --i;
i = low-1;
while ( ++i <= high ) {
p = x*V[i][k] + y*V[i][k+1];
if ( notlast ) { p += z*V[i][k+2]; V[i][k+2] -= p*r; }
......@@ -414,7 +420,7 @@ public class Eigen {
l = n; H[n][n] = 1; i = n;
while ( --i >= 0 ) {
w = H[i][i] - p;
r = 0; j = l; --j; while ( ++j <= n ) r += H[i][j] * H[j][n];
r = 0; j = l-1; while ( ++j <= n ) r += H[i][j] * H[j][n];
if ( e[i] < 0 ) { z = w; s = r; }
else {
l = i;
......@@ -425,60 +431,61 @@ public class Eigen {
}
// Overflow control
t = Math.abs(H[i][n]);
if ( EPS * t * t > 1 ) { j = i; --j; while ( ++j <= n ) H[j][n] /= t; }
if ( EPS * t * t > 1 ) { j = i-1; while ( ++j <= n ) H[j][n] /= t; }
}
}
}
else if ( q < 0 ) { // Complex vector
l = n-1;
// Last vector component imaginary so matrix is triangular
if ( Math.abs(H[n][n-1]) > Math.abs(H[n-1][n]) ) { H[n-1][n-1] = q / H[n][n-1]; H[n-1][n] = -(H[n][n] - p) / H[n][n-1]; }
else { cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q); H[n-1][n-1] = cdivr; H[n-1][n] = cdivi; } //################################ TO BE REIMPLEMENTED
H[n][n-1] = 0; H[n][n] = 1;
i = n-1;
while ( --i >= 0 ) {
ra = sa = 0;
j = l; --j; while ( ++j <= n ) { ra += H[i][j]*H[j][n-1]; sa += H[i][j]*H[j][n]; }
w = H[i][i] - p;
if ( e[i] < 0 ) { z = w; r = ra; s = sa; }
else {
l = i;
if ( e[i] == 0 ) {
cdiv(-ra,-sa,w,q); //################################ TO BE REIMPLEMENTED
H[i][n-1] = cdivr;
H[i][n] = cdivi;
}
else { // Solve complex equations
x = H[i][i+1];
y = H[i+1][i];
vr = square(d[i] - p) + square(e[i]) - q*q;
vi = (d[i] - p) * 2.0 * q;
if ( vr == 0.0 && vi == 0.0 ) vr = EPS * norm * (Math.abs(w) + Math.abs(q) + Math.abs(x) + Math.abs(y) + Math.abs(z));
cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi); //################################ TO BE REIMPLEMENTED
H[i][n-1] = cdivr;
H[i][n] = cdivi;
if ( Math.abs(x) > (Math.abs(z) + Math.abs(q)) ) { H[i+1][n-1] = (-ra - w*H[i][n-1] + q*H[i][n]) / x; H[i+1][n] = (-sa - w*H[i][n] - q*H[i][n-1]) / x; }
else {
cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q); //################################ TO BE REIMPLEMENTED
H[i+1][n-1] = cdivr;
H[i+1][n] = cdivi;
else
if ( q < 0 ) { // Complex vector
l = n-1;
// Last vector component imaginary so matrix is triangular
if ( Math.abs(H[n][n-1]) > Math.abs(H[n-1][n]) ) { H[n-1][n-1] = q / H[n][n-1]; H[n-1][n] = -(H[n][n] - p) / H[n][n-1]; }
else { cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q); H[n-1][n-1] = cdivr; H[n-1][n] = cdivi; } //################################ TO BE REIMPLEMENTED
H[n][n-1] = 0; H[n][n] = 1;
i = n-1;
while ( --i >= 0 ) {
ra = sa = 0;
j = l; --j; while ( ++j <= n ) { ra += H[i][j]*H[j][n-1]; sa += H[i][j]*H[j][n]; }
w = H[i][i] - p;
if ( e[i] < 0 ) { z = w; r = ra; s = sa; }
else {
l = i;
if ( e[i] == 0 ) {
cdiv(-ra,-sa,w,q); //################################ TO BE REIMPLEMENTED
H[i][n-1] = cdivr;
H[i][n] = cdivi;
}
else { // Solve complex equations
x = H[i][i+1];
y = H[i+1][i];
vr = square(d[i] - p) + square(e[i]) - q*q;
vi = (d[i] - p) * 2.0 * q;
if ( vr == 0.0 && vi == 0.0 ) vr = EPS * norm * (Math.abs(w) + Math.abs(q) + Math.abs(x) + Math.abs(y) + Math.abs(z));
cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi); //################################ TO BE REIMPLEMENTED
H[i][n-1] = cdivr;
H[i][n] = cdivi;
if ( Math.abs(x) > (Math.abs(z) + Math.abs(q)) ) { H[i+1][n-1] = (-ra - w*H[i][n-1] + q*H[i][n]) / x; H[i+1][n] = (-sa - w*H[i][n] - q*H[i][n-1]) / x; }
else {
cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q); //################################ TO BE REIMPLEMENTED
H[i+1][n-1] = cdivr;
H[i+1][n] = cdivi;
}
}
// Overflow control
t = Math.max(Math.abs(H[i][n-1]),Math.abs(H[i][n]));
if ( EPS * t * t > 1 ) { j = i; --j; while ( ++j <= n ) { H[j][n-1] /= t; H[j][n] /= t; } }
}
// Overflow control
t = Math.max(Math.abs(H[i][n-1]),Math.abs(H[i][n]));
if ( EPS * t * t > 1 ) { j = i; --j; while ( ++j <= n ) { H[j][n-1] /= t; H[j][n] /= t; } }
}
}
}
}
// Vectors of isolated roots
i = nn; while ( --i >= 0 ) if ( i < low || i > high ) { j = nn; while ( --j >= i ) V[i][j] = H[i][j]; }
for (int i: id0n) if ( i < low || i > high ) { j = i-1; while ( ++j < nn ) V[i][j] = H[i][j]; }
// Back transformation to get eigenvectors of original matrix
j = nn;
while ( --j >= low ) {
i = low; --i;
i = low-1;
while ( ++i <= high ) {
z = 0; k = ( j < high ) ? j+1 : high+1; while ( --k >= low ) z += V[i][k] * H[k][j];
z = 0; k = low-1; while ( ++k < ((j<high)?j:high) ) z += V[i][k] * H[k][j];
V[i][j] = z;
}
}
......@@ -489,12 +496,13 @@ public class Eigen {
//## Inversing eigenvector matrix V
private void inv() {
//## Initializing
nn = 2*this.n; H = new double[n][nn];
i = n; while ( --i >= 0 ) { H[i][n+i] = 1; j = n; while ( --j >= 0 ) H[i][j] = V[i][j]; }
nn = 2*this.n;
H = new double[n][nn];
for (int i: id0n ) { H[i][n+i] = 1; for (int j: id0n) H[i][j] = V[i][j]; }
//## Gauss-Jordan
i = n; while ( --i > 0 ) if ( H[i-1][0] < H[i][0] ) { j = nn; while ( --j >= 0 ) { r = H[i][j]; H[i][j] = H[i-1][j]; H[i-1][j] = r; } }
i = -1; while ( ++i < n ) { j = -1; while ( ++j < n ) if ( i != j ) { r = H[j][i]/H[i][i]; k = nn; while ( --k >= 0 ) H[j][k] -= r*H[i][k]; } }
i = -1; while ( ++i < n ) { r = H[i][i]; j = -1; while ( ++j < n ) inV[i][j] = H[i][n+j] / r; }
for (int i: id1n) if ( H[i-1][0] < H[i][0] ) { j = nn; while ( --j >= 0 ) { r = H[i][j]; H[i][j] = H[i-1][j]; H[i-1][j] = r; } }
for (int i: id0n) for (int j: id0n) if ( i != j ) { r = H[j][i]/H[i][i]; k = nn; while ( --k >= 0 ) H[j][k] -= r*H[i][k]; }
for (int i: id0n) { r = H[i][i]; for (int j: id0n) inV[i][j] = H[i][n+j] / r; }
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment