package math3d;
public class JacobiDouble { private double[][] matrix; private double[][] eigenmatrix; private double[] eigenvalues; private int numberOfRotationsNeeded; private int maxSweeps;
public JacobiDouble(double[][] matrix) { this(matrix, 50); }
public JacobiDouble(double[][] matrix, int maxSweeps) { this.matrix = matrix;
for(int i = 0; i < matrix.length; ++i) { for(int j = i + 1; j < matrix.length; ++j) { if (!this.isSmallComparedTo(Math.abs(matrix[i][j] - matrix[j][i]), matrix[i][j])) { throw new RuntimeException("Matrix is not symmetric!"); } } }
this.eigenmatrix = new double[matrix.length][matrix.length]; this.eigenvalues = new double[matrix.length]; this.maxSweeps = maxSweeps; this.perform(); }
public double[][] getEigenVectors() { return FastMatrixN.transpose(this.eigenmatrix); }
public double[][] getEigenMatrix() { return this.eigenmatrix; }
public double[] getEigenValues() { return this.eigenvalues; }
public int getNumberOfRotations() { return this.numberOfRotationsNeeded; }
private double offDiagonalSum() { double sum = 0.0;
for(int i = 0; i < this.matrix.length - 1; ++i) { for(int j = i + 1; j < this.matrix.length; ++j) { sum += Math.abs(this.matrix[i][j]); } }
return sum; }
private void rotate(int i, int j, int k, int l, double s, double tau) { double tmp1 = this.matrix[i][j]; double tmp2 = this.matrix[k][l]; this.matrix[i][j] = tmp1 - s * (tmp2 + tmp1 * tau); this.matrix[k][l] = tmp2 + s * (tmp1 - tmp2 * tau); }
private void rotateEigenMatrix(int i, int j, int k, int l, double s, double tau) { double tmp1 = this.eigenmatrix[i][j]; double tmp2 = this.eigenmatrix[k][l]; this.eigenmatrix[i][j] = tmp1 - s * (tmp2 + tmp1 * tau); this.eigenmatrix[k][l] = tmp2 + s * (tmp1 - tmp2 * tau); }
private boolean isSmallComparedTo(double value, double reference) { return Math.abs(reference) + value == Math.abs(reference); }
private void perform() { double[] b = new double[this.matrix.length]; double[] z = new double[this.matrix.length];
for(int i = 0; i < this.matrix.length; ++i) { for(int j = 0; j < this.matrix.length; ++j) { this.eigenmatrix[i][j] = 0.0; }
this.eigenmatrix[i][i] = 1.0; b[i] = this.eigenvalues[i] = this.matrix[i][i]; z[i] = 0.0; }
this.numberOfRotationsNeeded = 0;
for(int sweeps = 0; sweeps < this.maxSweeps; ++sweeps) { double sum = this.offDiagonalSum(); if (sum == 0.0) { return; }
double thresh = 0.0; if (sweeps < 3) { thresh = 0.2F * sum / (double)(this.matrix.length * this.matrix.length); }
for(int p = 0; p < this.matrix.length - 1; ++p) { for(int q = p + 1; q < this.matrix.length; ++q) { double tmp = 100.0 * Math.abs(this.matrix[p][q]); if (sweeps > 3 && this.isSmallComparedTo(tmp, this.eigenvalues[p]) && this.isSmallComparedTo(tmp, this.eigenvalues[q])) { this.matrix[p][q] = 0.0; } else if (Math.abs(this.matrix[p][q]) > thresh) { double diff = this.eigenvalues[q] - this.eigenvalues[p]; double t; if (this.isSmallComparedTo(tmp, diff)) { t = this.matrix[p][q] / diff; } else { double theta = 0.5 * diff / this.matrix[p][q]; t = 1.0 / (Math.abs(theta) + Math.sqrt(1.0 + theta * theta)); if (theta < 0.0) { t = -t; } }
double c = 1.0 / Math.sqrt(1.0 + t * t); double s = t * c; double tau = s / (1.0 + c); double h = t * this.matrix[p][q]; z[p] -= h; z[q] += h; this.eigenvalues[p] -= h; this.eigenvalues[q] += h; this.matrix[p][q] = 0.0;
for(int j = 0; j <= p - 1; ++j) { this.rotate(j, p, j, q, s, tau); }
for(int j = p + 1; j <= q - 1; ++j) { this.rotate(p, j, j, q, s, tau); }
for(int j = q + 1; j < this.matrix.length; ++j) { this.rotate(p, j, q, j, s, tau); }
for(int j = 0; j < this.matrix.length; ++j) { this.rotateEigenMatrix(j, p, j, q, s, tau); }
++this.numberOfRotationsNeeded; } } }
for(int p = 0; p < this.matrix.length; ++p) { b[p] += z[p]; this.eigenvalues[p] = b[p]; z[p] = 0.0; } } }
public static String toString(double[] doubleArray) { String result = "{";
for(int i = 0; i < doubleArray.length; ++i) { if (i > 0) { result = result + ","; }
result = result + doubleArray[i]; }
return result + "}"; }
public static String toString(double[][] double2Array) { String result = "{";
for(int i = 0; i < double2Array.length; ++i) { if (i > 0) { result = result + ","; }
result = result + toString(double2Array[i]); }
return result + "}"; }
public static double[] getColumn(double[][] matrix, int i) { double[] result = new double[matrix.length];
for(int j = 0; j < matrix.length; ++j) { result[j] = matrix[j][i]; }
return result; }
public static double[][] matMult(double[][] m1, double[][] m2) { int r = m1.length; int c = m2[0].length; double[][] result = new double[r][c];
for(int i = 0; i < r; ++i) { for(int j = 0; j < c; ++j) { result[i][j] = 0.0;
for(int k = 0; k < m2.length; ++k) { result[i][j] += m1[i][k] * m2[k][j]; } } }
return result; }
public static double[] vecMult(double[][] m, double[] v) { int r = m.length; double[] result = new double[r];
for(int i = 0; i < r; ++i) { result[i] = 0.0;
for(int k = 0; k < v.length; ++k) { result[i] += m[i][k] * v[k]; } }
return result; }
public static double[][] transpose(double[][] m) { int r = m.length; int c = m[0].length; double[][] result = new double[c][r];
for(int i = 0; i < r; ++i) { for(int j = 0; j < c; ++j) { result[j][i] = m[i][j]; } }
return result; }
public static void main(String[] args) { double[][] matrix = new double[][]{{1.0, 2.0}, {2.0, 1.0}}; JacobiDouble jacobi = new JacobiDouble(matrix); double[] eigenValuesVector = jacobi.getEigenValues(); double[][] eigenValues = new double[eigenValuesVector.length][eigenValuesVector.length];
for(int i = 0; i < eigenValuesVector.length; ++i) { eigenValues[i][i] = eigenValuesVector[i]; }
double[][] eigenVectors = jacobi.getEigenVectors(); double[][] result = matMult(eigenVectors, matMult(eigenValues, transpose(eigenVectors))); System.out.println("out: " + toString(result)); }}
