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.out.print("Too many iterations in tqli"); System.exit(0); } 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;k0) { 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;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