X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fmath%2FMatrix.java;h=b22bf4e1d27558fc09c6e59e64a45733f12d75cf;hb=187269da0130408d24901c88acb10c217f6284f9;hp=0f459788b5229bca917d5c4ddb5405d6477e48be;hpb=e5127e5f27f02c8f328539eaac68eabfbee1d135;p=jalview.git

diff --git a/src/jalview/math/Matrix.java b/src/jalview/math/Matrix.java
index 0f45978..b22bf4e 100755
--- a/src/jalview/math/Matrix.java
+++ b/src/jalview/math/Matrix.java
@@ -1,572 +1,1025 @@
-package jalview.math;
-
-import jalview.util.*;
-
-import java.io.*;
-
-public class Matrix {
-
-  /**
-   * SMJSPUBLIC
-   */
-  public double[][] value;
-  public int rows;
-  public int cols;
-  public double[] d;  // Diagonal
-  public double[] e;  // off diagonal
-
-  public Matrix(double[][] value, int rows, int cols) {
-    this.rows = rows;
-    this.cols = cols;
-    this.value = value;
-  }
-
-  public Matrix transpose() {
-    double[][] out = new double[cols][rows];
-
-    for (int i = 0; i < cols; i++) {
-      for (int j = 0; j < rows ; j++) {
-        out[i][j] = value[j][i];
-      }
-    }
-    return new Matrix(out,cols,rows);
-  }
-
-  public void print(PrintStream ps) {
-
-    for (int i = 0; i < rows; i++) {
-      for (int j = 0; j < cols; j++) {
-        Format.print(ps,"%8.2f",value[i][j]);
-      }
-      ps.println();
-    }
-  }
-
-
-  public Matrix  preMultiply(Matrix in) {
-    double[][] tmp = new double[in.rows][this.cols];
-
-    for (int i = 0; i < in.rows; i++) {
-      for (int j = 0; j < this.cols; j++ ) {
-        tmp[i][j] = 0.0;
-
-        for (int k = 0; k < in.cols; k++) {
-          tmp[i][j] += in.value[i][k]*this.value[k][j];
-        }
-
-      }
-    }
-
-    return new Matrix(tmp,in.rows,this.cols);
-  }
-
-  public double[] vectorPostMultiply(double[] in) {
-    double[] out = new double[in.length];
-    for (int i = 0; i < in.length; i++) {
-      out[i] = 0.0;
-      for (int k=0; k < in.length; k++) {
-        out[i] += value[i][k] * in[k];
-      }
-    }
-    return out;
-  }
-  public Matrix postMultiply(Matrix in) {
-
-    double[][] out = new double[this.rows][in.cols];
-    for (int i = 0; i < this.rows; i++) {
-      for (int j = 0; j < in.cols; j++ ) {
-
-        out[i][j] = 0.0;
-
-        for (int k = 0; k < rows; k++) {
-          out[i][j] = out[i][j] + value[i][k]*in.value[k][j];
-        }
-
-      }
-    }
-    return new Matrix(out,this.cols,in.rows);
-  }
-
-  public Matrix copy() {
-    double[][] newmat = new double[rows][cols];
-
-    for (int i = 0; i < rows; i++) {
-      for (int j = 0; j < cols; j++) {
-        newmat[i][j] = value[i][j];
-      }
-    }
-    return new Matrix(newmat,rows,cols);
-  }
-
-  public void tred() {
-    int n = rows;
-    int l;
-    int k;
-    int j;
-    int i;
-
-    double scale;
-    double hh;
-    double h;
-    double g;
-    double f;
-
-    this.d = new double[rows];
-    this.e = new double[rows];
-
-    for (i=n; i >= 2;i--) {
-      l=i-1;
-      h = 0.0;
-      scale = 0.0;
-
-      if (l > 1) {
-        for (k=1;k<=l;k++) {
-          scale += Math.abs(value[i-1][k-1]);
-        }
-        if (scale == 0.0) {
-          e[i-1] = value[i-1][l-1];
-        } else {
-          for (k=1; k <= l; k++) {
-            value[i-1][k-1] /= scale;
-            h += value[i-1][k-1]*value[i-1][k-1];
-          }
-          f = value[i-1][l-1];
-          if (f>0) {
-            g = -1.0*Math.sqrt(h);
-          } else {
-            g = Math.sqrt(h);
-          }
-          e[i-1] = scale*g;
-          h -= f*g;
-          value[i-1][l-1] = f-g;
-          f=0.0;
-          for (j=1; j <= l; j++) {
-            value[j-1][i-1] = value[i-1][j-1]/h;
-            g=0.0;
-            for (k= 1; k <= j; k++) {
-              g += value[j-1][k-1]*value[i-1][k-1];
-            }
-            for (k=j+1; k<=l;k++) {
-              g+= value[k-1][j-1]*value[i-1][k-1];
-            }
-            e[j-1] = g/h;
-            f+=e[j-1]*value[i-1][j-1];
-          }
-          hh=f/(h+h);
-          for (j=1;j<=l;j++) {
-            f=value[i-1][j-1];
-            g = e[j-1] - hh*f;
-            e[j-1] = g;
-            for (k=1;k<=j;k++) {
-              value[j-1][k-1] -= (f*e[k-1]+g*value[i-1][k-1]);
-            }
-          }
-        }
-      } else {
-        e[i-1] = value[i-1][l-1];
-      }
-      d[i-1] = h;
-    }
-    d[0] = 0.0;
-    e[0] = 0.0;
-    for (i=1;i<=n;i++) {
-      l=i-1;
-      if (d[i-1] != 0.0) {
-        for (j=1;j<=l;j++) {
-          g=0.0;
-          for (k=1;k<=l;k++) {
-            g+= value[i-1][k-1]*value[k-1][j-1];
-          }
-          for (k=1;k<=l;k++) {
-            value[k-1][j-1] -= g*value[k-1][i-1];
-          }
-        }
-      }
-      d[i-1] = value[i-1][i-1];
-      value[i-1][i-1] = 1.0;
-      for (j=1;j<=l;j++) {
-        value[j-1][i-1] = 0.0;
-        value[i-1][j-1] = 0.0;
-      }
-    }
-  }
-
-  public void tqli() {
-    int n = rows;
-
-    int m;
-    int l;
-    int iter;
-    int i;
-    int k;
-    double s;
-    double r;
-    double p;
-    ;
-    double g;
-    double f;
-    double dd;
-    double c;
-    double b;
-
-    for (i=2;i<=n;i++) {
-      e[i-2] = e[i-1];
-    }
-    e[n-1] = 0.0;
-    for (l=1;l<=n;l++) {
-      iter=0;
-      do {
-        for (m=l;m<=(n-1);m++) {
-          dd=Math.abs(d[m-1]) + Math.abs(d[m]);
-          if (Math.abs(e[m-1]) + dd == dd)
-            break;
-        }
-        if (m != l) {
-          iter++;
-          if (iter == 30) {
-            System.err.print("Too many iterations in tqli");
-            System.exit(0); // JBPNote - should this really be here ???
-          } else {
-            //	    System.out.println("Iteration " + iter);
-          }
-          g=(d[l]-d[l-1])/(2.0*e[l-1]);
-          r = Math.sqrt((g*g) + 1.0);
-          g=d[m-1]-d[l-1]+e[l-1]/(g + sign(r,g));
-          c = 1.0;
-          s = c;
-          p=0.0;
-          for (i=m-1;i>=l;i--) {
-            f = s*e[i-1];
-            b = c*e[i-1];
-            if (Math.abs(f) >= Math.abs(g)) {
-              c=g/f;
-              r = Math.sqrt((c*c)+1.0);
-              e[i] = f*r;
-              s = 1.0/r;
-              c *= s;
-            } else {
-              s=f/g;
-              r = Math.sqrt((s*s)+1.0);
-              e[i] = g*r;
-              c = 1.0/r;
-              s *= c;
-            }
-            g=d[i] -p;
-            r=(d[i-1]-g)*s + 2.0*c*b;
-            p=s*r;
-            d[i] = g + p;
-            g = c * r - b;
-            for (k=1; k <= n; k++) {
-              f=value[k-1][i];
-              value[k-1][i] = s*value[k-1][i-1] + c*f;
-              value[k-1][i-1] = c*value[k-1][i-1] - s*f;
-            }
-          }
-          d[l-1] = d[l-1] - p;
-          e[l-1] = g;
-          e[m-1] = 0.0;
-        }
-      } while ( m != l);
-    }
-  }
-  public void tred2() {
-    int n = rows;
-    int l;
-    int k;
-    int j;
-    int i;
-
-    double scale;
-    double hh;
-    double h;
-    double g;
-    double f;
-
-    this.d = new double[rows];
-    this.e = new double[rows];
-
-    for (i=n-1; i >= 1;i--) {
-      l=i-1;
-      h = 0.0;
-      scale = 0.0;
-
-      if (l > 0) {
-        for (k=0;k<l;k++) {
-          scale += Math.abs(value[i][k]);
-        }
-        if (scale == 0.0) {
-          e[i] = value[i][l];
-        } else {
-          for (k=0; k < l; k++) {
-            value[i][k] /= scale;
-            h += value[i][k]*value[i][k];
-          }
-          f = value[i][l];
-          if (f>0) {
-            g = -1.0*Math.sqrt(h);
-          } else {
-            g = Math.sqrt(h);
-          }
-          e[i] = scale*g;
-          h -= f*g;
-          value[i][l] = f-g;
-          f=0.0;
-          for (j=0; j < l; j++) {
-            value[j][i] = value[i][j]/h;
-            g=0.0;
-            for (k= 0; k < j; k++) {
-              g += value[j][k]*value[i][k];
-            }
-            for (k=j; k<l;k++) {
-              g+= value[k][j]*value[i][k];
-            }
-            e[j] = g/h;
-            f+=e[j]*value[i][j];
-          }
-          hh=f/(h+h);
-          for (j=0;j<l;j++) {
-            f=value[i][j];
-            g = e[j] - hh*f;
-            e[j] = g;
-            for (k=0;k<j;k++) {
-              value[j][k] -= (f*e[k]+g*value[i][k]);
-            }
-          }
-        }
-      } else {
-        e[i] = value[i][l];
-      }
-      d[i] = h;
-    }
-    d[0] = 0.0;
-    e[0] = 0.0;
-    for (i=0;i<n;i++) {
-      l=i-1;
-      if (d[i] != 0.0) {
-        for (j=0;j<l;j++) {
-          g=0.0;
-          for (k=0;k<l;k++) {
-            g+= value[i][k]*value[k][j];
-          }
-          for (k=0;k<l;k++) {
-            value[k][j] -= g*value[k][i];
-          }
-        }
-      }
-      d[i] = value[i][i];
-      value[i][i] = 1.0;
-      for (j=0;j<l;j++) {
-        value[j][i] = 0.0;
-        value[i][j] = 0.0;
-      }
-    }
-  }
-
-  public void tqli2() {
-    int n = rows;
-
-    int m;
-    int l;
-    int iter;
-    int i;
-    int k;
-    double s;
-    double r;
-    double p;
-    ;
-    double g;
-    double f;
-    double dd;
-    double c;
-    double b;
-
-    for (i=2;i<=n;i++) {
-      e[i-2] = e[i-1];
-    }
-    e[n-1] = 0.0;
-    for (l=1;l<=n;l++) {
-      iter=0;
-      do {
-        for (m=l;m<=(n-1);m++) {
-          dd=Math.abs(d[m-1]) + Math.abs(d[m]);
-          if (Math.abs(e[m-1]) + dd == dd)
-            break;
-        }
-        if (m != l) {
-          iter++;
-          if (iter == 30) {
-            System.err.print("Too many iterations in tqli");
-            System.exit(0); // JBPNote - same as above - not a graceful exit!
-          } else {
-            //	    System.out.println("Iteration " + iter);
-          }
-          g=(d[l]-d[l-1])/(2.0*e[l-1]);
-          r = Math.sqrt((g*g) + 1.0);
-          g=d[m-1]-d[l-1]+e[l-1]/(g + sign(r,g));
-          c = 1.0;
-          s = c;
-          p=0.0;
-          for (i=m-1;i>=l;i--) {
-            f = s*e[i-1];
-            b = c*e[i-1];
-            if (Math.abs(f) >= Math.abs(g)) {
-              c=g/f;
-              r = Math.sqrt((c*c)+1.0);
-              e[i] = f*r;
-              s = 1.0/r;
-              c *= s;
-            } else {
-              s=f/g;
-              r = Math.sqrt((s*s)+1.0);
-              e[i] = g*r;
-              c = 1.0/r;
-              s *= c;
-            }
-            g=d[i] -p;
-            r=(d[i-1]-g)*s + 2.0*c*b;
-            p=s*r;
-            d[i] = g + p;
-            g = c * r - b;
-            for (k=1; k <= n; k++) {
-              f=value[k-1][i];
-              value[k-1][i] = s*value[k-1][i-1] + c*f;
-              value[k-1][i-1] = c*value[k-1][i-1] - s*f;
-            }
-          }
-          d[l-1] = d[l-1] - p;
-          e[l-1] = g;
-          e[m-1] = 0.0;
-        }
-      } while ( m != l);
-    }
-  }
-
-  public double sign(double a, double b) {
-    if (b < 0) {
-      return -Math.abs(a);
-    } else {
-      return Math.abs(a);
-    }
-  }
-
-  public double[] getColumn(int n) {
-    double[] out  = new double[rows];
-    for (int i=0;i<rows;i++) {
-      out[i] = value[i][n];
-    }
-    return out;
-  }
-
-
-  public void printD(PrintStream ps) {
-
-    for (int j = 0; j < rows;j++) {
-      Format.print(ps,"%15.4e",d[j]);
-    }
-  }
-  public void printE(PrintStream ps) {
-
-    for (int j = 0; j < rows;j++) {
-      Format.print(ps,"%15.4e",e[j]);
-    }
-  }
-
-  public static void main(String[] args) {
-    int n = Integer.parseInt(args[0]);
-    double[][] in = new double[n][n];
-
-    for (int i = 0;i < n;i++) {
-      for (int j = 0; j < n; j++) {
-        in[i][j] = (double)Math.random();
-      }
-    }
-
-    Matrix origmat = new Matrix(in,n,n);
-    //    System.out.println(" --- Original matrix ---- ");
-    ///    origmat.print(System.out);
-    //System.out.println();
-
-    //System.out.println(" --- transpose matrix ---- ");
-    Matrix trans = origmat.transpose();
-    //trans.print(System.out);
-    //System.out.println();
-
-    //System.out.println(" --- OrigT * Orig ---- ");
-
-    Matrix symm = trans.postMultiply(origmat);
-    //symm.print(System.out);
-    //System.out.println();
-
-    // Copy the symmetric matrix for later
-    Matrix origsymm = symm.copy();
-
-
-    // This produces the tridiagonal transformation matrix
-    long tstart = System.currentTimeMillis();
-    symm.tred();
-    long tend = System.currentTimeMillis();
-    //System.out.println("Time take for tred = " + (tend-tstart) + "ms");
-    //System.out.println(" ---Tridiag transform matrix ---");
-    //symm.print(System.out);
-    //System.out.println();
-
-    //System.out.println(" --- D vector ---");
-    //symm.printD(System.out);
-    //System.out.println();
-    //System.out.println(" --- E vector ---");
-    //symm.printE(System.out);
-    //System.out.println();
-
-
-    // Now produce the diagonalization matrix
-    tstart = System.currentTimeMillis();
-    symm.tqli();
-    tend = System.currentTimeMillis();
-    //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
-
-    //System.out.println(" --- New diagonalization matrix ---");
-    //symm.print(System.out);
-    //System.out.println();
-
-    //System.out.println(" --- D vector ---");
-    //symm.printD(System.out);
-    //System.out.println();
-    //System.out.println(" --- E vector ---");
-    //symm.printE(System.out);
-    //System.out.println();
-
-    //System.out.println(" --- First eigenvector --- ");
-    //double[] eigenv = symm.getColumn(0);
-    //for (int i=0; i < eigenv.length;i++) {
-    //  Format.print(System.out,"%15.4f",eigenv[i]);
-    // }
-    //System.out.println();
-
-    //double[] neigenv = origsymm.vectorPostMultiply(eigenv);
-
-    //for (int i=0; i < neigenv.length;i++) {
-    //  Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);
-    //}
-
-    //System.out.println();
-  }
-
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+/*
+ * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
+ * Copyright (C) $$Year-Rel$$ The Jalview Authors
+ * 
+ * This file is part of Jalview.
+ * 
+ * Jalview is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License 
+ * as published by the Free Software Foundation, either version 3
+ * of the License, or (at your option) any later version.
+ *  
+ * Jalview is distributed in the hope that it will be useful, but 
+ * WITHOUT ANY WARRANTY; without even the implied warranty 
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
+ * PURPOSE.  See the GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
+ * The Jalview Authors are detailed in the 'AUTHORS' file.
+ */
+package jalview.math;
+
+import jalview.util.Format;
+import jalview.util.MessageManager;
+
+import java.io.PrintStream;
+import java.util.Arrays;
+
+/**
+ * A class to model rectangular matrices of double values and operations on them
+ */
+public class Matrix implements MatrixI
+{
+  /*
+   * maximum number of iterations for tqli
+   */
+  private static final int MAX_ITER = 45;
+  // fudge - add 15 iterations, just in case
+
+  /*
+   * the number of rows
+   */
+  final protected int rows;
+
+  /*
+   * the number of columns
+   */
+  final protected int cols;
+
+  /*
+   * the cell values in row-major order
+   */
+  private double[][] value;
+
+  protected double[] d; // Diagonal
+
+  protected double[] e; // off diagonal
+
+  /**
+   * Constructor given number of rows and columns
+   * 
+   * @param colCount
+   * @param rowCount
+   */
+  protected Matrix(int rowCount, int colCount)
+  {
+    rows = rowCount;
+    cols = colCount;
+  }
+
+  /**
+   * Creates a new Matrix object containing a copy of the supplied array values.
+   * For example
+   * 
+   * <pre>
+   *   new Matrix(new double[][] {{2, 3, 4}, {5, 6, 7})
+   * constructs
+   *   (2 3 4)
+   *   (5 6 7)
+   * </pre>
+   * 
+   * Note that ragged arrays (with not all rows, or columns, of the same
+   * length), are not supported by this class. They can be constructed, but
+   * results of operations on them are undefined and may throw exceptions.
+   * 
+   * @param values
+   *          the matrix values in row-major order
+   */
+  public Matrix(double[][] values)
+  {
+    this.rows = values.length;
+    this.cols = this.rows == 0 ? 0 : values[0].length;
+
+    /*
+     * make a copy of the values array, for immutability
+     */
+    this.value = new double[rows][];
+    int i = 0;
+    for (double[] row : values)
+    {
+      if (row != null)
+      {
+        value[i] = new double[row.length];
+        System.arraycopy(row, 0, value[i], 0, row.length);
+      }
+      i++;
+    }
+  }
+
+  @Override
+  public MatrixI transpose()
+  {
+    double[][] out = new double[cols][rows];
+
+    for (int i = 0; i < cols; i++)
+    {
+      for (int j = 0; j < rows; j++)
+      {
+        out[i][j] = value[j][i];
+      }
+    }
+
+    return new Matrix(out);
+  }
+
+  /**
+   * DOCUMENT ME!
+   * 
+   * @param ps
+   *          DOCUMENT ME!
+   * @param format
+   */
+  @Override
+  public void print(PrintStream ps, String format)
+  {
+    for (int i = 0; i < rows; i++)
+    {
+      for (int j = 0; j < cols; j++)
+      {
+        Format.print(ps, format, getValue(i, j));
+      }
+
+      ps.println();
+    }
+  }
+
+  @Override
+  public MatrixI preMultiply(MatrixI in)
+  {
+    if (in.width() != rows)
+    {
+      throw new IllegalArgumentException("Can't pre-multiply " + this.rows
+              + " rows by " + in.width() + " columns");
+    }
+    double[][] tmp = new double[in.height()][this.cols];
+
+    for (int i = 0; i < in.height(); i++)
+    {
+      for (int j = 0; j < this.cols; j++)
+      {
+        /*
+         * result[i][j] is the vector product of 
+         * in.row[i] and this.column[j]
+         */
+        for (int k = 0; k < in.width(); k++)
+        {
+          tmp[i][j] += (in.getValue(i, k) * this.value[k][j]);
+        }
+      }
+    }
+
+    return new Matrix(tmp);
+  }
+
+  /**
+   * 
+   * @param in
+   * 
+   * @return
+   */
+  public double[] vectorPostMultiply(double[] in)
+  {
+    double[] out = new double[in.length];
+
+    for (int i = 0; i < in.length; i++)
+    {
+      out[i] = 0.0;
+
+      for (int k = 0; k < in.length; k++)
+      {
+        out[i] += (value[i][k] * in[k]);
+      }
+    }
+
+    return out;
+  }
+
+  @Override
+  public MatrixI postMultiply(MatrixI in)
+  {
+    if (in.height() != this.cols)
+    {
+      throw new IllegalArgumentException("Can't post-multiply " + this.cols
+              + " columns by " + in.height() + " rows");
+    }
+    return in.preMultiply(this);
+  }
+
+  @Override
+  public MatrixI copy()
+  {
+    double[][] newmat = new double[rows][cols];
+
+    for (int i = 0; i < rows; i++)
+    {
+      System.arraycopy(value[i], 0, newmat[i], 0, value[i].length);
+    }
+
+    Matrix m = new Matrix(newmat);
+    if (this.d != null)
+    {
+      m.d = Arrays.copyOf(this.d, this.d.length);
+    }
+    if (this.e != null)
+    {
+      m.e = Arrays.copyOf(this.e, this.e.length);
+    }
+
+    return m;
+  }
+
+  /**
+   * DOCUMENT ME!
+   */
+  @Override
+  public void tred()
+  {
+    int n = rows;
+    int k;
+    int j;
+    int i;
+
+    double scale;
+    double hh;
+    double h;
+    double g;
+    double f;
+
+    this.d = new double[rows];
+    this.e = new double[rows];
+
+    for (i = n; i >= 2; i--)
+    {
+      final int l = i - 1;
+      h = 0.0;
+      scale = 0.0;
+
+      if (l > 1)
+      {
+        for (k = 1; k <= l; k++)
+        {
+          double v = Math.abs(getValue(i - 1, k - 1));
+          scale += v;
+        }
+
+        if (scale == 0.0)
+        {
+          e[i - 1] = getValue(i - 1, l - 1);
+        }
+        else
+        {
+          for (k = 1; k <= l; k++)
+          {
+            double v = divideValue(i - 1, k - 1, scale);
+            h += v * v;
+          }
+
+          f = getValue(i - 1, l - 1);
+
+          if (f > 0)
+          {
+            g = -1.0 * Math.sqrt(h);
+          }
+          else
+          {
+            g = Math.sqrt(h);
+          }
+
+          e[i - 1] = scale * g;
+          h -= (f * g);
+          setValue(i - 1, l - 1, f - g);
+          f = 0.0;
+
+          for (j = 1; j <= l; j++)
+          {
+            double val = getValue(i - 1, j - 1) / h;
+            setValue(j - 1, i - 1, val);
+            g = 0.0;
+
+            for (k = 1; k <= j; k++)
+            {
+              g += (getValue(j - 1, k - 1) * getValue(i - 1, k - 1));
+            }
+
+            for (k = j + 1; k <= l; k++)
+            {
+              g += (getValue(k - 1, j - 1) * getValue(i - 1, k - 1));
+            }
+
+            e[j - 1] = g / h;
+            f += (e[j - 1] * getValue(i - 1, j - 1));
+          }
+
+          hh = f / (h + h);
+
+          for (j = 1; j <= l; j++)
+          {
+            f = getValue(i - 1, j - 1);
+            g = e[j - 1] - (hh * f);
+            e[j - 1] = g;
+
+            for (k = 1; k <= j; k++)
+            {
+              double val = (f * e[k - 1]) + (g * getValue(i - 1, k - 1));
+              addValue(j - 1, k - 1, -val);
+            }
+          }
+        }
+      }
+      else
+      {
+        e[i - 1] = getValue(i - 1, l - 1);
+      }
+
+      d[i - 1] = h;
+    }
+
+    d[0] = 0.0;
+    e[0] = 0.0;
+
+    for (i = 1; i <= n; i++)
+    {
+      final int l = i - 1;
+
+      if (d[i - 1] != 0.0)
+      {
+        for (j = 1; j <= l; j++)
+        {
+          g = 0.0;
+
+          for (k = 1; k <= l; k++)
+          {
+            g += (getValue(i - 1, k - 1) * getValue(k - 1, j - 1));
+          }
+
+          for (k = 1; k <= l; k++)
+          {
+            addValue(k - 1, j - 1, -(g * getValue(k - 1, i - 1)));
+          }
+        }
+      }
+
+      d[i - 1] = getValue(i - 1, i - 1);
+      setValue(i - 1, i - 1, 1.0);
+
+      for (j = 1; j <= l; j++)
+      {
+        setValue(j - 1, i - 1, 0.0);
+        setValue(i - 1, j - 1, 0.0);
+      }
+    }
+  }
+
+  /**
+   * Adds f to the value at [i, j] and returns the new value
+   * 
+   * @param i
+   * @param j
+   * @param f
+   */
+  protected double addValue(int i, int j, double f)
+  {
+    double v = value[i][j] + f;
+    value[i][j] = v;
+    return v;
+  }
+
+  /**
+   * Divides the value at [i, j] by divisor and returns the new value. If d is
+   * zero, returns the unchanged value.
+   * 
+   * @param i
+   * @param j
+   * @param divisor
+   * @return
+   */
+  protected double divideValue(int i, int j, double divisor)
+  {
+    if (divisor == 0d)
+    {
+      return getValue(i, j);
+    }
+    double v = value[i][j];
+    v = v / divisor;
+    value[i][j] = v;
+    return v;
+  }
+
+  /**
+   * DOCUMENT ME!
+   */
+  @Override
+  public void tqli() throws Exception
+  {
+    int n = rows;
+
+    int m;
+    int l;
+    int iter;
+    int i;
+    int k;
+    double s;
+    double r;
+    double p;
+
+    double g;
+    double f;
+    double dd;
+    double c;
+    double b;
+
+    for (i = 2; i <= n; i++)
+    {
+      e[i - 2] = e[i - 1];
+    }
+
+    e[n - 1] = 0.0;
+
+    for (l = 1; l <= n; l++)
+    {
+      iter = 0;
+
+      do
+      {
+        for (m = l; m <= (n - 1); m++)
+        {
+          dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
+
+          if ((Math.abs(e[m - 1]) + dd) == dd)
+          {
+            break;
+          }
+        }
+
+        if (m != l)
+        {
+          iter++;
+
+          if (iter == MAX_ITER)
+          {
+            throw new Exception(MessageManager.formatMessage(
+                    "exception.matrix_too_many_iteration", new String[]
+                    { "tqli", Integer.valueOf(MAX_ITER).toString() }));
+          }
+          else
+          {
+            // System.out.println("Iteration " + iter);
+          }
+
+          g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
+          r = Math.sqrt((g * g) + 1.0);
+          g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
+          c = 1.0;
+          s = c;
+          p = 0.0;
+
+          for (i = m - 1; i >= l; i--)
+          {
+            f = s * e[i - 1];
+            b = c * e[i - 1];
+
+            if (Math.abs(f) >= Math.abs(g))
+            {
+              c = g / f;
+              r = Math.sqrt((c * c) + 1.0);
+              e[i] = f * r;
+              s = 1.0 / r;
+              c *= s;
+            }
+            else
+            {
+              s = f / g;
+              r = Math.sqrt((s * s) + 1.0);
+              e[i] = g * r;
+              c = 1.0 / r;
+              s *= c;
+            }
+
+            g = d[i] - p;
+            r = ((d[i - 1] - g) * s) + (2.0 * c * b);
+            p = s * r;
+            d[i] = g + p;
+            g = (c * r) - b;
+
+            for (k = 1; k <= n; k++)
+            {
+              f = getValue(k - 1, i);
+              setValue(k - 1, i, (s * getValue(k - 1, i - 1)) + (c * f));
+              setValue(k - 1, i - 1,
+                      (c * getValue(k - 1, i - 1)) - (s * f));
+            }
+          }
+
+          d[l - 1] = d[l - 1] - p;
+          e[l - 1] = g;
+          e[m - 1] = 0.0;
+        }
+      } while (m != l);
+    }
+  }
+
+  @Override
+  public double getValue(int i, int j)
+  {
+    return value[i][j];
+  }
+
+  @Override
+  public void setValue(int i, int j, double val)
+  {
+    value[i][j] = val;
+  }
+
+  /**
+   * DOCUMENT ME!
+   */
+  public void tred2()
+  {
+    int n = rows;
+    int l;
+    int k;
+    int j;
+    int i;
+
+    double scale;
+    double hh;
+    double h;
+    double g;
+    double f;
+
+    this.d = new double[rows];
+    this.e = new double[rows];
+
+    for (i = n - 1; i >= 1; i--)
+    {
+      l = i - 1;
+      h = 0.0;
+      scale = 0.0;
+
+      if (l > 0)
+      {
+        for (k = 0; k < l; k++)
+        {
+          scale += Math.abs(value[i][k]);
+        }
+
+        if (scale == 0.0)
+        {
+          e[i] = value[i][l];
+        }
+        else
+        {
+          for (k = 0; k < l; k++)
+          {
+            value[i][k] /= scale;
+            h += (value[i][k] * value[i][k]);
+          }
+
+          f = value[i][l];
+
+          if (f > 0)
+          {
+            g = -1.0 * Math.sqrt(h);
+          }
+          else
+          {
+            g = Math.sqrt(h);
+          }
+
+          e[i] = scale * g;
+          h -= (f * g);
+          value[i][l] = f - g;
+          f = 0.0;
+
+          for (j = 0; j < l; j++)
+          {
+            value[j][i] = value[i][j] / h;
+            g = 0.0;
+
+            for (k = 0; k < j; k++)
+            {
+              g += (value[j][k] * value[i][k]);
+            }
+
+            for (k = j; k < l; k++)
+            {
+              g += (value[k][j] * value[i][k]);
+            }
+
+            e[j] = g / h;
+            f += (e[j] * value[i][j]);
+          }
+
+          hh = f / (h + h);
+
+          for (j = 0; j < l; j++)
+          {
+            f = value[i][j];
+            g = e[j] - (hh * f);
+            e[j] = g;
+
+            for (k = 0; k < j; k++)
+            {
+              value[j][k] -= ((f * e[k]) + (g * value[i][k]));
+            }
+          }
+        }
+      }
+      else
+      {
+        e[i] = value[i][l];
+      }
+
+      d[i] = h;
+    }
+
+    d[0] = 0.0;
+    e[0] = 0.0;
+
+    for (i = 0; i < n; i++)
+    {
+      l = i - 1;
+
+      if (d[i] != 0.0)
+      {
+        for (j = 0; j < l; j++)
+        {
+          g = 0.0;
+
+          for (k = 0; k < l; k++)
+          {
+            g += (value[i][k] * value[k][j]);
+          }
+
+          for (k = 0; k < l; k++)
+          {
+            value[k][j] -= (g * value[k][i]);
+          }
+        }
+      }
+
+      d[i] = value[i][i];
+      value[i][i] = 1.0;
+
+      for (j = 0; j < l; j++)
+      {
+        value[j][i] = 0.0;
+        value[i][j] = 0.0;
+      }
+    }
+  }
+
+  /**
+   * DOCUMENT ME!
+   */
+  public void tqli2() throws Exception
+  {
+    int n = rows;
+
+    int m;
+    int l;
+    int iter;
+    int i;
+    int k;
+    double s;
+    double r;
+    double p;
+    ;
+
+    double g;
+    double f;
+    double dd;
+    double c;
+    double b;
+
+    for (i = 2; i <= n; i++)
+    {
+      e[i - 2] = e[i - 1];
+    }
+
+    e[n - 1] = 0.0;
+
+    for (l = 1; l <= n; l++)
+    {
+      iter = 0;
+
+      do
+      {
+        for (m = l; m <= (n - 1); m++)
+        {
+          dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
+
+          if ((Math.abs(e[m - 1]) + dd) == dd)
+          {
+            break;
+          }
+        }
+
+        if (m != l)
+        {
+          iter++;
+
+          if (iter == MAX_ITER)
+          {
+            throw new Exception(MessageManager.formatMessage(
+                    "exception.matrix_too_many_iteration", new String[]
+                    { "tqli2", Integer.valueOf(MAX_ITER).toString() }));
+          }
+          else
+          {
+            // System.out.println("Iteration " + iter);
+          }
+
+          g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
+          r = Math.sqrt((g * g) + 1.0);
+          g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
+          c = 1.0;
+          s = c;
+          p = 0.0;
+
+          for (i = m - 1; i >= l; i--)
+          {
+            f = s * e[i - 1];
+            b = c * e[i - 1];
+
+            if (Math.abs(f) >= Math.abs(g))
+            {
+              c = g / f;
+              r = Math.sqrt((c * c) + 1.0);
+              e[i] = f * r;
+              s = 1.0 / r;
+              c *= s;
+            }
+            else
+            {
+              s = f / g;
+              r = Math.sqrt((s * s) + 1.0);
+              e[i] = g * r;
+              c = 1.0 / r;
+              s *= c;
+            }
+
+            g = d[i] - p;
+            r = ((d[i - 1] - g) * s) + (2.0 * c * b);
+            p = s * r;
+            d[i] = g + p;
+            g = (c * r) - b;
+
+            for (k = 1; k <= n; k++)
+            {
+              f = value[k - 1][i];
+              value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
+              value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
+            }
+          }
+
+          d[l - 1] = d[l - 1] - p;
+          e[l - 1] = g;
+          e[m - 1] = 0.0;
+        }
+      } while (m != l);
+    }
+  }
+
+  /**
+   * Answers the first argument with the sign of the second argument
+   * 
+   * @param a
+   * @param b
+   * 
+   * @return
+   */
+  static double sign(double a, double b)
+  {
+    if (b < 0)
+    {
+      return -Math.abs(a);
+    }
+    else
+    {
+      return Math.abs(a);
+    }
+  }
+
+  /**
+   * Returns an array containing the values in the specified column
+   * 
+   * @param col
+   * 
+   * @return
+   */
+  public double[] getColumn(int col)
+  {
+    double[] out = new double[rows];
+
+    for (int i = 0; i < rows; i++)
+    {
+      out[i] = value[i][col];
+    }
+
+    return out;
+  }
+
+  /**
+   * DOCUMENT ME!
+   * 
+   * @param ps
+   *          DOCUMENT ME!
+   * @param format
+   */
+  @Override
+  public void printD(PrintStream ps, String format)
+  {
+    for (int j = 0; j < rows; j++)
+    {
+      Format.print(ps, format, d[j]);
+    }
+  }
+
+  /**
+   * DOCUMENT ME!
+   * 
+   * @param ps
+   *          DOCUMENT ME!
+   * @param format
+   *          TODO
+   */
+  @Override
+  public void printE(PrintStream ps, String format)
+  {
+    for (int j = 0; j < rows; j++)
+    {
+      Format.print(ps, format, e[j]);
+    }
+  }
+
+  @Override
+  public double[] getD()
+  {
+    return d;
+  }
+
+  @Override
+  public double[] getE()
+  {
+    return e;
+  }
+
+  @Override
+  public int height()
+  {
+    return rows;
+  }
+
+  @Override
+  public int width()
+  {
+    return cols;
+  }
+
+  @Override
+  public double[] getRow(int i)
+  {
+    double[] row = new double[cols];
+    System.arraycopy(value[i], 0, row, 0, cols);
+    return row;
+  }
+
+  /**
+   * Returns a length 2 array of {minValue, maxValue} of all values in the
+   * matrix. Returns null if the matrix is null or empty.
+   * 
+   * @return
+   */
+  double[] findMinMax()
+  {
+    if (value == null)
+    {
+      return null;
+    }
+    double min = Double.MAX_VALUE;
+    double max = -Double.MAX_VALUE;
+    boolean empty = true;
+    for (double[] row : value)
+    {
+      if (row != null)
+      {
+        for (double x : row)
+        {
+          empty = false;
+          if (x > max)
+          {
+            max = x;
+          }
+          if (x < min)
+          {
+            min = x;
+          }
+        }
+      }
+    }
+    return empty ? null : new double[] { min, max };
+  }
+
+  /**
+   * {@inheritDoc}
+   */
+  @Override
+  public void reverseRange(boolean maxToZero)
+  {
+    if (value == null)
+    {
+      return;
+    }
+    double[] minMax = findMinMax();
+    if (minMax == null)
+    {
+      return; // empty matrix
+    }
+    double subtractFrom = maxToZero ? minMax[1] : minMax[0] + minMax[1];
+
+    for (double[] row : value)
+    {
+      if (row != null)
+      {
+        int j = 0;
+        for (double x : row)
+        {
+          row[j] = subtractFrom - x;
+          j++;
+        }
+      }
+    }
+  }
+
+  /**
+   * Multiplies every entry in the matrix by the given value.
+   * 
+   * @param
+   */
+  @Override
+  public void multiply(double by)
+  {
+    for (double[] row : value)
+    {
+      if (row != null)
+      {
+        for (int i = 0; i < row.length; i++)
+        {
+          row[i] *= by;
+        }
+      }
+    }
+  }
+
+  @Override
+  public void setD(double[] v)
+  {
+    d = v;
+  }
+
+  @Override
+  public void setE(double[] v)
+  {
+    e = v;
+  }
+
+  public double getTotal()
+  {
+    double d = 0d;
+    for (int i = 0; i < this.height(); i++)
+    {
+      for (int j = 0; j < this.width(); j++)
+      {
+        d += value[i][j];
+      }
+    }
+    return d;
+  }
+
+  /**
+   * {@inheritDoc}
+   */
+  @Override
+  public boolean equals(MatrixI m2, double delta)
+  {
+    if (m2 == null || this.height() != m2.height()
+            || this.width() != m2.width())
+    {
+      return false;
+    }
+    for (int i = 0; i < this.height(); i++)
+    {
+      for (int j = 0; j < this.width(); j++)
+      {
+        double diff = this.getValue(i, j) - m2.getValue(i, j);
+        if (Math.abs(diff) > delta)
+        {
+          return false;
+        }
+      }
+    }
+    return true;
+  }
+}