src/jalview/math/Matrix.java

   1 package jalview.math;\r
   2 \r
   3 import jalview.util.*;\r
   4 \r
   5 import java.io.*;\r
   6 \r
   7 public class Matrix {\r
   8 \r
   9   /**\r
  10    * SMJSPUBLIC\r
  11    */\r
  12   public double[][] value;\r
  13   public int rows;\r
  14   public int cols;\r
  15   public double[] d;  // Diagonal\r
  16   public double[] e;  // off diagonal\r
  17 \r
  18   public Matrix(double[][] value, int rows, int cols) {\r
  19     this.rows = rows;\r
  20     this.cols = cols;\r
  21     this.value = value;\r
  22   }\r
  23 \r
  24   public Matrix transpose() {\r
  25     double[][] out = new double[cols][rows];\r
  26 \r
  27     for (int i = 0; i < cols; i++) {\r
  28       for (int j = 0; j < rows ; j++) {\r
  29         out[i][j] = value[j][i];\r
  30       }\r
  31     }\r
  32     return new Matrix(out,cols,rows);\r
  33   }\r
  34 \r
  35   public void print(PrintStream ps) {\r
  36 \r
  37     for (int i = 0; i < rows; i++) {\r
  38       for (int j = 0; j < cols; j++) {\r
  39         Format.print(ps,"%8.2f",value[i][j]);\r
  40       }\r
  41       ps.println();\r
  42     }\r
  43   }\r
  44 \r
  45 \r
  46   public Matrix  preMultiply(Matrix in) {\r
  47     double[][] tmp = new double[in.rows][this.cols];\r
  48 \r
  49     for (int i = 0; i < in.rows; i++) {\r
  50       for (int j = 0; j < this.cols; j++ ) {\r
  51         tmp[i][j] = 0.0;\r
  52 \r
  53         for (int k = 0; k < in.cols; k++) {\r
  54           tmp[i][j] += in.value[i][k]*this.value[k][j];\r
  55         }\r
  56 \r
  57       }\r
  58     }\r
  59 \r
  60     return new Matrix(tmp,in.rows,this.cols);\r
  61   }\r
  62 \r
  63   public double[] vectorPostMultiply(double[] in) {\r
  64     double[] out = new double[in.length];\r
  65     for (int i = 0; i < in.length; i++) {\r
  66       out[i] = 0.0;\r
  67       for (int k=0; k < in.length; k++) {\r
  68         out[i] += value[i][k] * in[k];\r
  69       }\r
  70     }\r
  71     return out;\r
  72   }\r
  73   public Matrix postMultiply(Matrix in) {\r
  74 \r
  75     double[][] out = new double[this.rows][in.cols];\r
  76     for (int i = 0; i < this.rows; i++) {\r
  77       for (int j = 0; j < in.cols; j++ ) {\r
  78 \r
  79         out[i][j] = 0.0;\r
  80 \r
  81         for (int k = 0; k < rows; k++) {\r
  82           out[i][j] = out[i][j] + value[i][k]*in.value[k][j];\r
  83         }\r
  84 \r
  85       }\r
  86     }\r
  87     return new Matrix(out,this.cols,in.rows);\r
  88   }\r
  89 \r
  90   public Matrix copy() {\r
  91     double[][] newmat = new double[rows][cols];\r
  92 \r
  93     for (int i = 0; i < rows; i++) {\r
  94       for (int j = 0; j < cols; j++) {\r
  95         newmat[i][j] = value[i][j];\r
  96       }\r
  97     }\r
  98     return new Matrix(newmat,rows,cols);\r
  99   }\r
 100 \r
 101   public void tred() {\r
 102     int n = rows;\r
 103     int l;\r
 104     int k;\r
 105     int j;\r
 106     int i;\r
 107 \r
 108     double scale;\r
 109     double hh;\r
 110     double h;\r
 111     double g;\r
 112     double f;\r
 113 \r
 114     this.d = new double[rows];\r
 115     this.e = new double[rows];\r
 116 \r
 117     for (i=n; i >= 2;i--) {\r
 118       l=i-1;\r
 119       h = 0.0;\r
 120       scale = 0.0;\r
 121 \r
 122       if (l > 1) {\r
 123         for (k=1;k<=l;k++) {\r
 124           scale += Math.abs(value[i-1][k-1]);\r
 125         }\r
 126         if (scale == 0.0) {\r
 127           e[i-1] = value[i-1][l-1];\r
 128         } else {\r
 129           for (k=1; k <= l; k++) {\r
 130             value[i-1][k-1] /= scale;\r
 131             h += value[i-1][k-1]*value[i-1][k-1];\r
 132           }\r
 133           f = value[i-1][l-1];\r
 134           if (f>0) {\r
 135             g = -1.0*Math.sqrt(h);\r
 136           } else {\r
 137             g = Math.sqrt(h);\r
 138           }\r
 139           e[i-1] = scale*g;\r
 140           h -= f*g;\r
 141           value[i-1][l-1] = f-g;\r
 142           f=0.0;\r
 143           for (j=1; j <= l; j++) {\r
 144             value[j-1][i-1] = value[i-1][j-1]/h;\r
 145             g=0.0;\r
 146             for (k= 1; k <= j; k++) {\r
 147               g += value[j-1][k-1]*value[i-1][k-1];\r
 148             }\r
 149             for (k=j+1; k<=l;k++) {\r
 150               g+= value[k-1][j-1]*value[i-1][k-1];\r
 151             }\r
 152             e[j-1] = g/h;\r
 153             f+=e[j-1]*value[i-1][j-1];\r
 154           }\r
 155           hh=f/(h+h);\r
 156           for (j=1;j<=l;j++) {\r
 157             f=value[i-1][j-1];\r
 158             g = e[j-1] - hh*f;\r
 159             e[j-1] = g;\r
 160             for (k=1;k<=j;k++) {\r
 161               value[j-1][k-1] -= (f*e[k-1]+g*value[i-1][k-1]);\r
 162             }\r
 163           }\r
 164         }\r
 165       } else {\r
 166         e[i-1] = value[i-1][l-1];\r
 167       }\r
 168       d[i-1] = h;\r
 169     }\r
 170     d[0] = 0.0;\r
 171     e[0] = 0.0;\r
 172     for (i=1;i<=n;i++) {\r
 173       l=i-1;\r
 174       if (d[i-1] != 0.0) {\r
 175         for (j=1;j<=l;j++) {\r
 176           g=0.0;\r
 177           for (k=1;k<=l;k++) {\r
 178             g+= value[i-1][k-1]*value[k-1][j-1];\r
 179           }\r
 180           for (k=1;k<=l;k++) {\r
 181             value[k-1][j-1] -= g*value[k-1][i-1];\r
 182           }\r
 183         }\r
 184       }\r
 185       d[i-1] = value[i-1][i-1];\r
 186       value[i-1][i-1] = 1.0;\r
 187       for (j=1;j<=l;j++) {\r
 188         value[j-1][i-1] = 0.0;\r
 189         value[i-1][j-1] = 0.0;\r
 190       }\r
 191     }\r
 192   }\r
 193 \r
 194   public void tqli() {\r
 195     int n = rows;\r
 196 \r
 197     int m;\r
 198     int l;\r
 199     int iter;\r
 200     int i;\r
 201     int k;\r
 202     double s;\r
 203     double r;\r
 204     double p;\r
 205     ;\r
 206     double g;\r
 207     double f;\r
 208     double dd;\r
 209     double c;\r
 210     double b;\r
 211 \r
 212     for (i=2;i<=n;i++) {\r
 213       e[i-2] = e[i-1];\r
 214     }\r
 215     e[n-1] = 0.0;\r
 216     for (l=1;l<=n;l++) {\r
 217       iter=0;\r
 218       do {\r
 219         for (m=l;m<=(n-1);m++) {\r
 220           dd=Math.abs(d[m-1]) + Math.abs(d[m]);\r
 221           if (Math.abs(e[m-1]) + dd == dd)\r
 222             break;\r
 223         }\r
 224         if (m != l) {\r
 225           iter++;\r
 226           if (iter == 30) {\r
 227             System.err.print("Too many iterations in tqli");\r
 228             System.exit(0); // JBPNote - should this really be here ???\r
 229           } else {\r
 230             //      System.out.println("Iteration " + iter);\r
 231           }\r
 232           g=(d[l]-d[l-1])/(2.0*e[l-1]);\r
 233           r = Math.sqrt((g*g) + 1.0);\r
 234           g=d[m-1]-d[l-1]+e[l-1]/(g + sign(r,g));\r
 235           c = 1.0;\r
 236           s = c;\r
 237           p=0.0;\r
 238           for (i=m-1;i>=l;i--) {\r
 239             f = s*e[i-1];\r
 240             b = c*e[i-1];\r
 241             if (Math.abs(f) >= Math.abs(g)) {\r
 242               c=g/f;\r
 243               r = Math.sqrt((c*c)+1.0);\r
 244               e[i] = f*r;\r
 245               s = 1.0/r;\r
 246               c *= s;\r
 247             } else {\r
 248               s=f/g;\r
 249               r = Math.sqrt((s*s)+1.0);\r
 250               e[i] = g*r;\r
 251               c = 1.0/r;\r
 252               s *= c;\r
 253             }\r
 254             g=d[i] -p;\r
 255             r=(d[i-1]-g)*s + 2.0*c*b;\r
 256             p=s*r;\r
 257             d[i] = g + p;\r
 258             g = c * r - b;\r
 259             for (k=1; k <= n; k++) {\r
 260               f=value[k-1][i];\r
 261               value[k-1][i] = s*value[k-1][i-1] + c*f;\r
 262               value[k-1][i-1] = c*value[k-1][i-1] - s*f;\r
 263             }\r
 264           }\r
 265           d[l-1] = d[l-1] - p;\r
 266           e[l-1] = g;\r
 267           e[m-1] = 0.0;\r
 268         }\r
 269       } while ( m != l);\r
 270     }\r
 271   }\r
 272   public void tred2() {\r
 273     int n = rows;\r
 274     int l;\r
 275     int k;\r
 276     int j;\r
 277     int i;\r
 278 \r
 279     double scale;\r
 280     double hh;\r
 281     double h;\r
 282     double g;\r
 283     double f;\r
 284 \r
 285     this.d = new double[rows];\r
 286     this.e = new double[rows];\r
 287 \r
 288     for (i=n-1; i >= 1;i--) {\r
 289       l=i-1;\r
 290       h = 0.0;\r
 291       scale = 0.0;\r
 292 \r
 293       if (l > 0) {\r
 294         for (k=0;k<l;k++) {\r
 295           scale += Math.abs(value[i][k]);\r
 296         }\r
 297         if (scale == 0.0) {\r
 298           e[i] = value[i][l];\r
 299         } else {\r
 300           for (k=0; k < l; k++) {\r
 301             value[i][k] /= scale;\r
 302             h += value[i][k]*value[i][k];\r
 303           }\r
 304           f = value[i][l];\r
 305           if (f>0) {\r
 306             g = -1.0*Math.sqrt(h);\r
 307           } else {\r
 308             g = Math.sqrt(h);\r
 309           }\r
 310           e[i] = scale*g;\r
 311           h -= f*g;\r
 312           value[i][l] = f-g;\r
 313           f=0.0;\r
 314           for (j=0; j < l; j++) {\r
 315             value[j][i] = value[i][j]/h;\r
 316             g=0.0;\r
 317             for (k= 0; k < j; k++) {\r
 318               g += value[j][k]*value[i][k];\r
 319             }\r
 320             for (k=j; k<l;k++) {\r
 321               g+= value[k][j]*value[i][k];\r
 322             }\r
 323             e[j] = g/h;\r
 324             f+=e[j]*value[i][j];\r
 325           }\r
 326           hh=f/(h+h);\r
 327           for (j=0;j<l;j++) {\r
 328             f=value[i][j];\r
 329             g = e[j] - hh*f;\r
 330             e[j] = g;\r
 331             for (k=0;k<j;k++) {\r
 332               value[j][k] -= (f*e[k]+g*value[i][k]);\r
 333             }\r
 334           }\r
 335         }\r
 336       } else {\r
 337         e[i] = value[i][l];\r
 338       }\r
 339       d[i] = h;\r
 340     }\r
 341     d[0] = 0.0;\r
 342     e[0] = 0.0;\r
 343     for (i=0;i<n;i++) {\r
 344       l=i-1;\r
 345       if (d[i] != 0.0) {\r
 346         for (j=0;j<l;j++) {\r
 347           g=0.0;\r
 348           for (k=0;k<l;k++) {\r
 349             g+= value[i][k]*value[k][j];\r
 350           }\r
 351           for (k=0;k<l;k++) {\r
 352             value[k][j] -= g*value[k][i];\r
 353           }\r
 354         }\r
 355       }\r
 356       d[i] = value[i][i];\r
 357       value[i][i] = 1.0;\r
 358       for (j=0;j<l;j++) {\r
 359         value[j][i] = 0.0;\r
 360         value[i][j] = 0.0;\r
 361       }\r
 362     }\r
 363   }\r
 364 \r
 365   public void tqli2() {\r
 366     int n = rows;\r
 367 \r
 368     int m;\r
 369     int l;\r
 370     int iter;\r
 371     int i;\r
 372     int k;\r
 373     double s;\r
 374     double r;\r
 375     double p;\r
 376     ;\r
 377     double g;\r
 378     double f;\r
 379     double dd;\r
 380     double c;\r
 381     double b;\r
 382 \r
 383     for (i=2;i<=n;i++) {\r
 384       e[i-2] = e[i-1];\r
 385     }\r
 386     e[n-1] = 0.0;\r
 387     for (l=1;l<=n;l++) {\r
 388       iter=0;\r
 389       do {\r
 390         for (m=l;m<=(n-1);m++) {\r
 391           dd=Math.abs(d[m-1]) + Math.abs(d[m]);\r
 392           if (Math.abs(e[m-1]) + dd == dd)\r
 393             break;\r
 394         }\r
 395         if (m != l) {\r
 396           iter++;\r
 397           if (iter == 30) {\r
 398             System.err.print("Too many iterations in tqli");\r
 399             System.exit(0); // JBPNote - same as above - not a graceful exit!\r
 400           } else {\r
 401             //      System.out.println("Iteration " + iter);\r
 402           }\r
 403           g=(d[l]-d[l-1])/(2.0*e[l-1]);\r
 404           r = Math.sqrt((g*g) + 1.0);\r
 405           g=d[m-1]-d[l-1]+e[l-1]/(g + sign(r,g));\r
 406           c = 1.0;\r
 407           s = c;\r
 408           p=0.0;\r
 409           for (i=m-1;i>=l;i--) {\r
 410             f = s*e[i-1];\r
 411             b = c*e[i-1];\r
 412             if (Math.abs(f) >= Math.abs(g)) {\r
 413               c=g/f;\r
 414               r = Math.sqrt((c*c)+1.0);\r
 415               e[i] = f*r;\r
 416               s = 1.0/r;\r
 417               c *= s;\r
 418             } else {\r
 419               s=f/g;\r
 420               r = Math.sqrt((s*s)+1.0);\r
 421               e[i] = g*r;\r
 422               c = 1.0/r;\r
 423               s *= c;\r
 424             }\r
 425             g=d[i] -p;\r
 426             r=(d[i-1]-g)*s + 2.0*c*b;\r
 427             p=s*r;\r
 428             d[i] = g + p;\r
 429             g = c * r - b;\r
 430             for (k=1; k <= n; k++) {\r
 431               f=value[k-1][i];\r
 432               value[k-1][i] = s*value[k-1][i-1] + c*f;\r
 433               value[k-1][i-1] = c*value[k-1][i-1] - s*f;\r
 434             }\r
 435           }\r
 436           d[l-1] = d[l-1] - p;\r
 437           e[l-1] = g;\r
 438           e[m-1] = 0.0;\r
 439         }\r
 440       } while ( m != l);\r
 441     }\r
 442   }\r
 443 \r
 444   public double sign(double a, double b) {\r
 445     if (b < 0) {\r
 446       return -Math.abs(a);\r
 447     } else {\r
 448       return Math.abs(a);\r
 449     }\r
 450   }\r
 451 \r
 452   public double[] getColumn(int n) {\r
 453     double[] out  = new double[rows];\r
 454     for (int i=0;i<rows;i++) {\r
 455       out[i] = value[i][n];\r
 456     }\r
 457     return out;\r
 458   }\r
 459 \r
 460 \r
 461   public void printD(PrintStream ps) {\r
 462 \r
 463     for (int j = 0; j < rows;j++) {\r
 464       Format.print(ps,"%15.4e",d[j]);\r
 465     }\r
 466   }\r
 467   public void printE(PrintStream ps) {\r
 468 \r
 469     for (int j = 0; j < rows;j++) {\r
 470       Format.print(ps,"%15.4e",e[j]);\r
 471     }\r
 472   }\r
 473 \r
 474   public static void main(String[] args) {\r
 475     int n = Integer.parseInt(args[0]);\r
 476     double[][] in = new double[n][n];\r
 477 \r
 478     for (int i = 0;i < n;i++) {\r
 479       for (int j = 0; j < n; j++) {\r
 480         in[i][j] = (double)Math.random();\r
 481       }\r
 482     }\r
 483 \r
 484     Matrix origmat = new Matrix(in,n,n);\r
 485     //    System.out.println(" --- Original matrix ---- ");\r
 486     ///    origmat.print(System.out);\r
 487     //System.out.println();\r
 488 \r
 489     //System.out.println(" --- transpose matrix ---- ");\r
 490     Matrix trans = origmat.transpose();\r
 491     //trans.print(System.out);\r
 492     //System.out.println();\r
 493 \r
 494     //System.out.println(" --- OrigT * Orig ---- ");\r
 495 \r
 496     Matrix symm = trans.postMultiply(origmat);\r
 497     //symm.print(System.out);\r
 498     //System.out.println();\r
 499 \r
 500     // Copy the symmetric matrix for later\r
 501     Matrix origsymm = symm.copy();\r
 502 \r
 503 \r
 504     // This produces the tridiagonal transformation matrix\r
 505     long tstart = System.currentTimeMillis();\r
 506     symm.tred();\r
 507     long tend = System.currentTimeMillis();\r
 508     //System.out.println("Time take for tred = " + (tend-tstart) + "ms");\r
 509     //System.out.println(" ---Tridiag transform matrix ---");\r
 510     //symm.print(System.out);\r
 511     //System.out.println();\r
 512 \r
 513     //System.out.println(" --- D vector ---");\r
 514     //symm.printD(System.out);\r
 515     //System.out.println();\r
 516     //System.out.println(" --- E vector ---");\r
 517     //symm.printE(System.out);\r
 518     //System.out.println();\r
 519 \r
 520 \r
 521     // Now produce the diagonalization matrix\r
 522     tstart = System.currentTimeMillis();\r
 523     symm.tqli();\r
 524     tend = System.currentTimeMillis();\r
 525     //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");\r
 526 \r
 527     //System.out.println(" --- New diagonalization matrix ---");\r
 528     //symm.print(System.out);\r
 529     //System.out.println();\r
 530 \r
 531     //System.out.println(" --- D vector ---");\r
 532     //symm.printD(System.out);\r
 533     //System.out.println();\r
 534     //System.out.println(" --- E vector ---");\r
 535     //symm.printE(System.out);\r
 536     //System.out.println();\r
 537 \r
 538     //System.out.println(" --- First eigenvector --- ");\r
 539     //double[] eigenv = symm.getColumn(0);\r
 540     //for (int i=0; i < eigenv.length;i++) {\r
 541     //  Format.print(System.out,"%15.4f",eigenv[i]);\r
 542     // }\r
 543     //System.out.println();\r
 544 \r
 545     //double[] neigenv = origsymm.vectorPostMultiply(eigenv);\r
 546 \r
 547     //for (int i=0; i < neigenv.length;i++) {\r
 548     //  Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);\r
 549     //}\r
 550 \r
 551     //System.out.println();\r
 552   }\r
 553 \r
 554 }\r
 555 \r
 556 \r
 557 \r
 558 \r
 559 \r
 560 \r
 561 \r
 562 \r
 563 \r
 564 \r
 565 \r
 566 \r
 567 \r
 568 \r
 569 \r
 570 \r
 571 \r
 572 \r