src/jalview/math/Matrix.java

   1 /*\r
   2 * Jalview - A Sequence Alignment Editor and Viewer\r
   3 * Copyright (C) 2006 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle\r
   4 *\r
   5 * This program is free software; you can redistribute it and/or\r
   6 * modify it under the terms of the GNU General Public License\r
   7 * as published by the Free Software Foundation; either version 2\r
   8 * of the License, or (at your option) any later version.\r
   9 *\r
  10 * This program is distributed in the hope that it will be useful,\r
  11 * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  13 * GNU General Public License for more details.\r
  14 *\r
  15 * You should have received a copy of the GNU General Public License\r
  16 * along with this program; if not, write to the Free Software\r
  17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA\r
  18 */\r
  19 package jalview.math;\r
  20 \r
  21 import jalview.util.*;\r
  22 \r
  23 import java.io.*;\r
  24 \r
  25 \r
  26 /**\r
  27  * DOCUMENT ME!\r
  28  *\r
  29  * @author $author$\r
  30  * @version $Revision$\r
  31  */\r
  32 public class Matrix\r
  33 {\r
  34     /**\r
  35      * SMJSPUBLIC\r
  36      */\r
  37     public double[][] value;\r
  38 \r
  39     /** DOCUMENT ME!! */\r
  40     public int rows;\r
  41 \r
  42     /** DOCUMENT ME!! */\r
  43     public int cols;\r
  44 \r
  45     /** DOCUMENT ME!! */\r
  46     public double[] d; // Diagonal\r
  47 \r
  48     /** DOCUMENT ME!! */\r
  49     public double[] e; // off diagonal\r
  50 \r
  51     /**\r
  52      * Creates a new Matrix object.\r
  53      *\r
  54      * @param value DOCUMENT ME!\r
  55      * @param rows DOCUMENT ME!\r
  56      * @param cols DOCUMENT ME!\r
  57      */\r
  58     public Matrix(double[][] value, int rows, int cols)\r
  59     {\r
  60         this.rows = rows;\r
  61         this.cols = cols;\r
  62         this.value = value;\r
  63     }\r
  64 \r
  65     /**\r
  66      * DOCUMENT ME!\r
  67      *\r
  68      * @return DOCUMENT ME!\r
  69      */\r
  70     public Matrix transpose()\r
  71     {\r
  72         double[][] out = new double[cols][rows];\r
  73 \r
  74         for (int i = 0; i < cols; i++)\r
  75         {\r
  76             for (int j = 0; j < rows; j++)\r
  77             {\r
  78                 out[i][j] = value[j][i];\r
  79             }\r
  80         }\r
  81 \r
  82         return new Matrix(out, cols, rows);\r
  83     }\r
  84 \r
  85     /**\r
  86      * DOCUMENT ME!\r
  87      *\r
  88      * @param ps DOCUMENT ME!\r
  89      */\r
  90     public void print(PrintStream ps)\r
  91     {\r
  92         for (int i = 0; i < rows; i++)\r
  93         {\r
  94             for (int j = 0; j < cols; j++)\r
  95             {\r
  96                 Format.print(ps, "%8.2f", value[i][j]);\r
  97             }\r
  98 \r
  99             ps.println();\r
 100         }\r
 101     }\r
 102 \r
 103     /**\r
 104      * DOCUMENT ME!\r
 105      *\r
 106      * @param in DOCUMENT ME!\r
 107      *\r
 108      * @return DOCUMENT ME!\r
 109      */\r
 110     public Matrix preMultiply(Matrix in)\r
 111     {\r
 112         double[][] tmp = new double[in.rows][this.cols];\r
 113 \r
 114         for (int i = 0; i < in.rows; i++)\r
 115         {\r
 116             for (int j = 0; j < this.cols; j++)\r
 117             {\r
 118                 tmp[i][j] = 0.0;\r
 119 \r
 120                 for (int k = 0; k < in.cols; k++)\r
 121                 {\r
 122                     tmp[i][j] += (in.value[i][k] * this.value[k][j]);\r
 123                 }\r
 124             }\r
 125         }\r
 126 \r
 127         return new Matrix(tmp, in.rows, this.cols);\r
 128     }\r
 129 \r
 130     /**\r
 131      * DOCUMENT ME!\r
 132      *\r
 133      * @param in DOCUMENT ME!\r
 134      *\r
 135      * @return DOCUMENT ME!\r
 136      */\r
 137     public double[] vectorPostMultiply(double[] in)\r
 138     {\r
 139         double[] out = new double[in.length];\r
 140 \r
 141         for (int i = 0; i < in.length; i++)\r
 142         {\r
 143             out[i] = 0.0;\r
 144 \r
 145             for (int k = 0; k < in.length; k++)\r
 146             {\r
 147                 out[i] += (value[i][k] * in[k]);\r
 148             }\r
 149         }\r
 150 \r
 151         return out;\r
 152     }\r
 153 \r
 154     /**\r
 155      * DOCUMENT ME!\r
 156      *\r
 157      * @param in DOCUMENT ME!\r
 158      *\r
 159      * @return DOCUMENT ME!\r
 160      */\r
 161     public Matrix postMultiply(Matrix in)\r
 162     {\r
 163         double[][] out = new double[this.rows][in.cols];\r
 164 \r
 165         for (int i = 0; i < this.rows; i++)\r
 166         {\r
 167             for (int j = 0; j < in.cols; j++)\r
 168             {\r
 169                 out[i][j] = 0.0;\r
 170 \r
 171                 for (int k = 0; k < rows; k++)\r
 172                 {\r
 173                     out[i][j] = out[i][j] + (value[i][k] * in.value[k][j]);\r
 174                 }\r
 175             }\r
 176         }\r
 177 \r
 178         return new Matrix(out, this.cols, in.rows);\r
 179     }\r
 180 \r
 181     /**\r
 182      * DOCUMENT ME!\r
 183      *\r
 184      * @return DOCUMENT ME!\r
 185      */\r
 186     public Matrix copy()\r
 187     {\r
 188         double[][] newmat = new double[rows][cols];\r
 189 \r
 190         for (int i = 0; i < rows; i++)\r
 191         {\r
 192             for (int j = 0; j < cols; j++)\r
 193             {\r
 194                 newmat[i][j] = value[i][j];\r
 195             }\r
 196         }\r
 197 \r
 198         return new Matrix(newmat, rows, cols);\r
 199     }\r
 200 \r
 201     /**\r
 202      * DOCUMENT ME!\r
 203      */\r
 204     public void tred()\r
 205     {\r
 206         int n = rows;\r
 207         int l;\r
 208         int k;\r
 209         int j;\r
 210         int i;\r
 211 \r
 212         double scale;\r
 213         double hh;\r
 214         double h;\r
 215         double g;\r
 216         double f;\r
 217 \r
 218         this.d = new double[rows];\r
 219         this.e = new double[rows];\r
 220 \r
 221         for (i = n; i >= 2; i--)\r
 222         {\r
 223             l = i - 1;\r
 224             h = 0.0;\r
 225             scale = 0.0;\r
 226 \r
 227             if (l > 1)\r
 228             {\r
 229                 for (k = 1; k <= l; k++)\r
 230                 {\r
 231                     scale += Math.abs(value[i - 1][k - 1]);\r
 232                 }\r
 233 \r
 234                 if (scale == 0.0)\r
 235                 {\r
 236                     e[i - 1] = value[i - 1][l - 1];\r
 237                 }\r
 238                 else\r
 239                 {\r
 240                     for (k = 1; k <= l; k++)\r
 241                     {\r
 242                         value[i - 1][k - 1] /= scale;\r
 243                         h += (value[i - 1][k - 1] * value[i - 1][k - 1]);\r
 244                     }\r
 245 \r
 246                     f = value[i - 1][l - 1];\r
 247 \r
 248                     if (f > 0)\r
 249                     {\r
 250                         g = -1.0 * Math.sqrt(h);\r
 251                     }\r
 252                     else\r
 253                     {\r
 254                         g = Math.sqrt(h);\r
 255                     }\r
 256 \r
 257                     e[i - 1] = scale * g;\r
 258                     h -= (f * g);\r
 259                     value[i - 1][l - 1] = f - g;\r
 260                     f = 0.0;\r
 261 \r
 262                     for (j = 1; j <= l; j++)\r
 263                     {\r
 264                         value[j - 1][i - 1] = value[i - 1][j - 1] / h;\r
 265                         g = 0.0;\r
 266 \r
 267                         for (k = 1; k <= j; k++)\r
 268                         {\r
 269                             g += (value[j - 1][k - 1] * value[i - 1][k - 1]);\r
 270                         }\r
 271 \r
 272                         for (k = j + 1; k <= l; k++)\r
 273                         {\r
 274                             g += (value[k - 1][j - 1] * value[i - 1][k - 1]);\r
 275                         }\r
 276 \r
 277                         e[j - 1] = g / h;\r
 278                         f += (e[j - 1] * value[i - 1][j - 1]);\r
 279                     }\r
 280 \r
 281                     hh = f / (h + h);\r
 282 \r
 283                     for (j = 1; j <= l; j++)\r
 284                     {\r
 285                         f = value[i - 1][j - 1];\r
 286                         g = e[j - 1] - (hh * f);\r
 287                         e[j - 1] = g;\r
 288 \r
 289                         for (k = 1; k <= j; k++)\r
 290                         {\r
 291                             value[j - 1][k - 1] -= ((f * e[k - 1]) +\r
 292                             (g * value[i - 1][k - 1]));\r
 293                         }\r
 294                     }\r
 295                 }\r
 296             }\r
 297             else\r
 298             {\r
 299                 e[i - 1] = value[i - 1][l - 1];\r
 300             }\r
 301 \r
 302             d[i - 1] = h;\r
 303         }\r
 304 \r
 305         d[0] = 0.0;\r
 306         e[0] = 0.0;\r
 307 \r
 308         for (i = 1; i <= n; i++)\r
 309         {\r
 310             l = i - 1;\r
 311 \r
 312             if (d[i - 1] != 0.0)\r
 313             {\r
 314                 for (j = 1; j <= l; j++)\r
 315                 {\r
 316                     g = 0.0;\r
 317 \r
 318                     for (k = 1; k <= l; k++)\r
 319                     {\r
 320                         g += (value[i - 1][k - 1] * value[k - 1][j - 1]);\r
 321                     }\r
 322 \r
 323                     for (k = 1; k <= l; k++)\r
 324                     {\r
 325                         value[k - 1][j - 1] -= (g * value[k - 1][i - 1]);\r
 326                     }\r
 327                 }\r
 328             }\r
 329 \r
 330             d[i - 1] = value[i - 1][i - 1];\r
 331             value[i - 1][i - 1] = 1.0;\r
 332 \r
 333             for (j = 1; j <= l; j++)\r
 334             {\r
 335                 value[j - 1][i - 1] = 0.0;\r
 336                 value[i - 1][j - 1] = 0.0;\r
 337             }\r
 338         }\r
 339     }\r
 340 \r
 341     /**\r
 342      * DOCUMENT ME!\r
 343      */\r
 344     public void tqli()\r
 345     {\r
 346         int n = rows;\r
 347 \r
 348         int m;\r
 349         int l;\r
 350         int iter;\r
 351         int i;\r
 352         int k;\r
 353         double s;\r
 354         double r;\r
 355         double p;\r
 356         ;\r
 357 \r
 358         double g;\r
 359         double f;\r
 360         double dd;\r
 361         double c;\r
 362         double b;\r
 363 \r
 364         for (i = 2; i <= n; i++)\r
 365         {\r
 366             e[i - 2] = e[i - 1];\r
 367         }\r
 368 \r
 369         e[n - 1] = 0.0;\r
 370 \r
 371         for (l = 1; l <= n; l++)\r
 372         {\r
 373             iter = 0;\r
 374 \r
 375             do\r
 376             {\r
 377                 for (m = l; m <= (n - 1); m++)\r
 378                 {\r
 379                     dd = Math.abs(d[m - 1]) + Math.abs(d[m]);\r
 380 \r
 381                     if ((Math.abs(e[m - 1]) + dd) == dd)\r
 382                     {\r
 383                         break;\r
 384                     }\r
 385                 }\r
 386 \r
 387                 if (m != l)\r
 388                 {\r
 389                     iter++;\r
 390 \r
 391                     if (iter == 30)\r
 392                     {\r
 393                         System.err.print("Too many iterations in tqli");\r
 394                         System.exit(0); // JBPNote - should this really be here ???\r
 395                     }\r
 396                     else\r
 397                     {\r
 398                         //          System.out.println("Iteration " + iter);\r
 399                     }\r
 400 \r
 401                     g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);\r
 402                     r = Math.sqrt((g * g) + 1.0);\r
 403                     g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));\r
 404                     c = 1.0;\r
 405                     s = c;\r
 406                     p = 0.0;\r
 407 \r
 408                     for (i = m - 1; i >= l; i--)\r
 409                     {\r
 410                         f = s * e[i - 1];\r
 411                         b = c * e[i - 1];\r
 412 \r
 413                         if (Math.abs(f) >= Math.abs(g))\r
 414                         {\r
 415                             c = g / f;\r
 416                             r = Math.sqrt((c * c) + 1.0);\r
 417                             e[i] = f * r;\r
 418                             s = 1.0 / r;\r
 419                             c *= s;\r
 420                         }\r
 421                         else\r
 422                         {\r
 423                             s = f / g;\r
 424                             r = Math.sqrt((s * s) + 1.0);\r
 425                             e[i] = g * r;\r
 426                             c = 1.0 / r;\r
 427                             s *= c;\r
 428                         }\r
 429 \r
 430                         g = d[i] - p;\r
 431                         r = ((d[i - 1] - g) * s) + (2.0 * c * b);\r
 432                         p = s * r;\r
 433                         d[i] = g + p;\r
 434                         g = (c * r) - b;\r
 435 \r
 436                         for (k = 1; k <= n; k++)\r
 437                         {\r
 438                             f = value[k - 1][i];\r
 439                             value[k - 1][i] = (s * value[k - 1][i - 1]) +\r
 440                                 (c * f);\r
 441                             value[k - 1][i - 1] = (c * value[k - 1][i - 1]) -\r
 442                                 (s * f);\r
 443                         }\r
 444                     }\r
 445 \r
 446                     d[l - 1] = d[l - 1] - p;\r
 447                     e[l - 1] = g;\r
 448                     e[m - 1] = 0.0;\r
 449                 }\r
 450             }\r
 451             while (m != l);\r
 452         }\r
 453     }\r
 454 \r
 455     /**\r
 456      * DOCUMENT ME!\r
 457      */\r
 458     public void tred2()\r
 459     {\r
 460         int n = rows;\r
 461         int l;\r
 462         int k;\r
 463         int j;\r
 464         int i;\r
 465 \r
 466         double scale;\r
 467         double hh;\r
 468         double h;\r
 469         double g;\r
 470         double f;\r
 471 \r
 472         this.d = new double[rows];\r
 473         this.e = new double[rows];\r
 474 \r
 475         for (i = n - 1; i >= 1; i--)\r
 476         {\r
 477             l = i - 1;\r
 478             h = 0.0;\r
 479             scale = 0.0;\r
 480 \r
 481             if (l > 0)\r
 482             {\r
 483                 for (k = 0; k < l; k++)\r
 484                 {\r
 485                     scale += Math.abs(value[i][k]);\r
 486                 }\r
 487 \r
 488                 if (scale == 0.0)\r
 489                 {\r
 490                     e[i] = value[i][l];\r
 491                 }\r
 492                 else\r
 493                 {\r
 494                     for (k = 0; k < l; k++)\r
 495                     {\r
 496                         value[i][k] /= scale;\r
 497                         h += (value[i][k] * value[i][k]);\r
 498                     }\r
 499 \r
 500                     f = value[i][l];\r
 501 \r
 502                     if (f > 0)\r
 503                     {\r
 504                         g = -1.0 * Math.sqrt(h);\r
 505                     }\r
 506                     else\r
 507                     {\r
 508                         g = Math.sqrt(h);\r
 509                     }\r
 510 \r
 511                     e[i] = scale * g;\r
 512                     h -= (f * g);\r
 513                     value[i][l] = f - g;\r
 514                     f = 0.0;\r
 515 \r
 516                     for (j = 0; j < l; j++)\r
 517                     {\r
 518                         value[j][i] = value[i][j] / h;\r
 519                         g = 0.0;\r
 520 \r
 521                         for (k = 0; k < j; k++)\r
 522                         {\r
 523                             g += (value[j][k] * value[i][k]);\r
 524                         }\r
 525 \r
 526                         for (k = j; k < l; k++)\r
 527                         {\r
 528                             g += (value[k][j] * value[i][k]);\r
 529                         }\r
 530 \r
 531                         e[j] = g / h;\r
 532                         f += (e[j] * value[i][j]);\r
 533                     }\r
 534 \r
 535                     hh = f / (h + h);\r
 536 \r
 537                     for (j = 0; j < l; j++)\r
 538                     {\r
 539                         f = value[i][j];\r
 540                         g = e[j] - (hh * f);\r
 541                         e[j] = g;\r
 542 \r
 543                         for (k = 0; k < j; k++)\r
 544                         {\r
 545                             value[j][k] -= ((f * e[k]) + (g * value[i][k]));\r
 546                         }\r
 547                     }\r
 548                 }\r
 549             }\r
 550             else\r
 551             {\r
 552                 e[i] = value[i][l];\r
 553             }\r
 554 \r
 555             d[i] = h;\r
 556         }\r
 557 \r
 558         d[0] = 0.0;\r
 559         e[0] = 0.0;\r
 560 \r
 561         for (i = 0; i < n; i++)\r
 562         {\r
 563             l = i - 1;\r
 564 \r
 565             if (d[i] != 0.0)\r
 566             {\r
 567                 for (j = 0; j < l; j++)\r
 568                 {\r
 569                     g = 0.0;\r
 570 \r
 571                     for (k = 0; k < l; k++)\r
 572                     {\r
 573                         g += (value[i][k] * value[k][j]);\r
 574                     }\r
 575 \r
 576                     for (k = 0; k < l; k++)\r
 577                     {\r
 578                         value[k][j] -= (g * value[k][i]);\r
 579                     }\r
 580                 }\r
 581             }\r
 582 \r
 583             d[i] = value[i][i];\r
 584             value[i][i] = 1.0;\r
 585 \r
 586             for (j = 0; j < l; j++)\r
 587             {\r
 588                 value[j][i] = 0.0;\r
 589                 value[i][j] = 0.0;\r
 590             }\r
 591         }\r
 592     }\r
 593 \r
 594     /**\r
 595      * DOCUMENT ME!\r
 596      */\r
 597     public void tqli2()\r
 598     {\r
 599         int n = rows;\r
 600 \r
 601         int m;\r
 602         int l;\r
 603         int iter;\r
 604         int i;\r
 605         int k;\r
 606         double s;\r
 607         double r;\r
 608         double p;\r
 609         ;\r
 610 \r
 611         double g;\r
 612         double f;\r
 613         double dd;\r
 614         double c;\r
 615         double b;\r
 616 \r
 617         for (i = 2; i <= n; i++)\r
 618         {\r
 619             e[i - 2] = e[i - 1];\r
 620         }\r
 621 \r
 622         e[n - 1] = 0.0;\r
 623 \r
 624         for (l = 1; l <= n; l++)\r
 625         {\r
 626             iter = 0;\r
 627 \r
 628             do\r
 629             {\r
 630                 for (m = l; m <= (n - 1); m++)\r
 631                 {\r
 632                     dd = Math.abs(d[m - 1]) + Math.abs(d[m]);\r
 633 \r
 634                     if ((Math.abs(e[m - 1]) + dd) == dd)\r
 635                     {\r
 636                         break;\r
 637                     }\r
 638                 }\r
 639 \r
 640                 if (m != l)\r
 641                 {\r
 642                     iter++;\r
 643 \r
 644                     if (iter == 30)\r
 645                     {\r
 646                         System.err.print("Too many iterations in tqli");\r
 647                         System.exit(0); // JBPNote - same as above - not a graceful exit!\r
 648                     }\r
 649                     else\r
 650                     {\r
 651                         //          System.out.println("Iteration " + iter);\r
 652                     }\r
 653 \r
 654                     g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);\r
 655                     r = Math.sqrt((g * g) + 1.0);\r
 656                     g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));\r
 657                     c = 1.0;\r
 658                     s = c;\r
 659                     p = 0.0;\r
 660 \r
 661                     for (i = m - 1; i >= l; i--)\r
 662                     {\r
 663                         f = s * e[i - 1];\r
 664                         b = c * e[i - 1];\r
 665 \r
 666                         if (Math.abs(f) >= Math.abs(g))\r
 667                         {\r
 668                             c = g / f;\r
 669                             r = Math.sqrt((c * c) + 1.0);\r
 670                             e[i] = f * r;\r
 671                             s = 1.0 / r;\r
 672                             c *= s;\r
 673                         }\r
 674                         else\r
 675                         {\r
 676                             s = f / g;\r
 677                             r = Math.sqrt((s * s) + 1.0);\r
 678                             e[i] = g * r;\r
 679                             c = 1.0 / r;\r
 680                             s *= c;\r
 681                         }\r
 682 \r
 683                         g = d[i] - p;\r
 684                         r = ((d[i - 1] - g) * s) + (2.0 * c * b);\r
 685                         p = s * r;\r
 686                         d[i] = g + p;\r
 687                         g = (c * r) - b;\r
 688 \r
 689                         for (k = 1; k <= n; k++)\r
 690                         {\r
 691                             f = value[k - 1][i];\r
 692                             value[k - 1][i] = (s * value[k - 1][i - 1]) +\r
 693                                 (c * f);\r
 694                             value[k - 1][i - 1] = (c * value[k - 1][i - 1]) -\r
 695                                 (s * f);\r
 696                         }\r
 697                     }\r
 698 \r
 699                     d[l - 1] = d[l - 1] - p;\r
 700                     e[l - 1] = g;\r
 701                     e[m - 1] = 0.0;\r
 702                 }\r
 703             }\r
 704             while (m != l);\r
 705         }\r
 706     }\r
 707 \r
 708     /**\r
 709      * DOCUMENT ME!\r
 710      *\r
 711      * @param a DOCUMENT ME!\r
 712      * @param b DOCUMENT ME!\r
 713      *\r
 714      * @return DOCUMENT ME!\r
 715      */\r
 716     public double sign(double a, double b)\r
 717     {\r
 718         if (b < 0)\r
 719         {\r
 720             return -Math.abs(a);\r
 721         }\r
 722         else\r
 723         {\r
 724             return Math.abs(a);\r
 725         }\r
 726     }\r
 727 \r
 728     /**\r
 729      * DOCUMENT ME!\r
 730      *\r
 731      * @param n DOCUMENT ME!\r
 732      *\r
 733      * @return DOCUMENT ME!\r
 734      */\r
 735     public double[] getColumn(int n)\r
 736     {\r
 737         double[] out = new double[rows];\r
 738 \r
 739         for (int i = 0; i < rows; i++)\r
 740         {\r
 741             out[i] = value[i][n];\r
 742         }\r
 743 \r
 744         return out;\r
 745     }\r
 746 \r
 747     /**\r
 748      * DOCUMENT ME!\r
 749      *\r
 750      * @param ps DOCUMENT ME!\r
 751      */\r
 752     public void printD(PrintStream ps)\r
 753     {\r
 754         for (int j = 0; j < rows; j++)\r
 755         {\r
 756             Format.print(ps, "%15.4e", d[j]);\r
 757         }\r
 758     }\r
 759 \r
 760     /**\r
 761      * DOCUMENT ME!\r
 762      *\r
 763      * @param ps DOCUMENT ME!\r
 764      */\r
 765     public void printE(PrintStream ps)\r
 766     {\r
 767         for (int j = 0; j < rows; j++)\r
 768         {\r
 769             Format.print(ps, "%15.4e", e[j]);\r
 770         }\r
 771     }\r
 772 \r
 773     /**\r
 774      * DOCUMENT ME!\r
 775      *\r
 776      * @param args DOCUMENT ME!\r
 777      */\r
 778     public static void main(String[] args)\r
 779     {\r
 780         int n = Integer.parseInt(args[0]);\r
 781         double[][] in = new double[n][n];\r
 782 \r
 783         for (int i = 0; i < n; i++)\r
 784         {\r
 785             for (int j = 0; j < n; j++)\r
 786             {\r
 787                 in[i][j] = (double) Math.random();\r
 788             }\r
 789         }\r
 790 \r
 791         Matrix origmat = new Matrix(in, n, n);\r
 792 \r
 793         //    System.out.println(" --- Original matrix ---- ");\r
 794         ///    origmat.print(System.out);\r
 795         //System.out.println();\r
 796         //System.out.println(" --- transpose matrix ---- ");\r
 797         Matrix trans = origmat.transpose();\r
 798 \r
 799         //trans.print(System.out);\r
 800         //System.out.println();\r
 801         //System.out.println(" --- OrigT * Orig ---- ");\r
 802         Matrix symm = trans.postMultiply(origmat);\r
 803 \r
 804         //symm.print(System.out);\r
 805         //System.out.println();\r
 806         // Copy the symmetric matrix for later\r
 807         //Matrix origsymm = symm.copy();\r
 808 \r
 809         // This produces the tridiagonal transformation matrix\r
 810         //long tstart = System.currentTimeMillis();\r
 811         symm.tred();\r
 812 \r
 813         //long tend = System.currentTimeMillis();\r
 814 \r
 815         //System.out.println("Time take for tred = " + (tend-tstart) + "ms");\r
 816         //System.out.println(" ---Tridiag transform matrix ---");\r
 817         //symm.print(System.out);\r
 818         //System.out.println();\r
 819         //System.out.println(" --- D vector ---");\r
 820         //symm.printD(System.out);\r
 821         //System.out.println();\r
 822         //System.out.println(" --- E vector ---");\r
 823         //symm.printE(System.out);\r
 824         //System.out.println();\r
 825         // Now produce the diagonalization matrix\r
 826         //tstart = System.currentTimeMillis();\r
 827         symm.tqli();\r
 828         //tend = System.currentTimeMillis();\r
 829 \r
 830         //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");\r
 831         //System.out.println(" --- New diagonalization matrix ---");\r
 832         //symm.print(System.out);\r
 833         //System.out.println();\r
 834         //System.out.println(" --- D vector ---");\r
 835         //symm.printD(System.out);\r
 836         //System.out.println();\r
 837         //System.out.println(" --- E vector ---");\r
 838         //symm.printE(System.out);\r
 839         //System.out.println();\r
 840         //System.out.println(" --- First eigenvector --- ");\r
 841         //double[] eigenv = symm.getColumn(0);\r
 842         //for (int i=0; i < eigenv.length;i++) {\r
 843         //  Format.print(System.out,"%15.4f",eigenv[i]);\r
 844         // }\r
 845         //System.out.println();\r
 846         //double[] neigenv = origsymm.vectorPostMultiply(eigenv);\r
 847         //for (int i=0; i < neigenv.length;i++) {\r
 848         //  Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);\r
 849         //}\r
 850         //System.out.println();\r
 851     }\r
 852 }\r