src/javajs/util/Matrix.java

   1 package javajs.util;\r
   2 \r
   3 /**\r
   4  * \r
   5  * streamlined and refined for Jmol by Bob Hanson\r
   6  * \r
   7  * from http://math.nist.gov/javanumerics/jama/\r
   8  * \r
   9  * Jama = Java Matrix class.\r
  10  * \r
  11  * @author The MathWorks, Inc. and the National Institute of Standards and\r
  12  *         Technology.\r
  13  * @version 5 August 1998\r
  14  */\r
  15 \r
  16 public class Matrix implements Cloneable {\r
  17 \r
  18   public double[][] a;\r
  19   protected int m, n;\r
  20 \r
  21   /**\r
  22    * Construct a matrix quickly without checking arguments.\r
  23    * \r
  24    * @param a\r
  25    *        Two-dimensional array of doubles or null\r
  26    * @param m\r
  27    *        Number of rows.\r
  28    * @param n\r
  29    *        Number of colums.\r
  30    */\r
  31 \r
  32   public Matrix(double[][] a, int m, int n) {\r
  33     this.a = (a == null ? new double[m][n] : a);\r
  34     this.m = m;\r
  35     this.n = n;\r
  36   }\r
  37 \r
  38   /**\r
  39    * Get row dimension.\r
  40    * \r
  41    * @return m, the number of rows.\r
  42    */\r
  43 \r
  44   public int getRowDimension() {\r
  45     return m;\r
  46   }\r
  47 \r
  48   /**\r
  49    * Get column dimension.\r
  50    * \r
  51    * @return n, the number of columns.\r
  52    */\r
  53 \r
  54   public int getColumnDimension() {\r
  55     return n;\r
  56   }\r
  57 \r
  58   /**\r
  59    * Access the internal two-dimensional array.\r
  60    * \r
  61    * @return Pointer to the two-dimensional array of matrix elements.\r
  62    */\r
  63 \r
  64   public double[][] getArray() {\r
  65     return a;\r
  66   }\r
  67 \r
  68   /**\r
  69    * Copy the internal two-dimensional array.\r
  70    * \r
  71    * @return Two-dimensional array copy of matrix elements.\r
  72    */\r
  73 \r
  74   public double[][] getArrayCopy() {\r
  75     double[][] x = new double[m][n];\r
  76     for (int i = m; --i >= 0;)\r
  77       for (int j = n; --j >= 0;)\r
  78         x[i][j] = a[i][j];\r
  79     return x;\r
  80   }\r
  81 \r
  82   /**\r
  83    * Make a deep copy of a matrix\r
  84    * \r
  85    * @return copy\r
  86    */\r
  87 \r
  88   public Matrix copy() {\r
  89     Matrix x = new Matrix(null, m, n);\r
  90     double[][] c = x.a;\r
  91     for (int i = m; --i >= 0;)\r
  92       for (int j = n; --j >= 0;)\r
  93         c[i][j] = a[i][j];\r
  94     return x;\r
  95   }\r
  96 \r
  97   /**\r
  98    * Clone the Matrix object.\r
  99    */\r
 100 \r
 101   @Override\r
 102   public Object clone() {\r
 103     return copy();\r
 104   }\r
 105 \r
 106   /**\r
 107    * Get a submatrix.\r
 108    * \r
 109    * @param i0\r
 110    *        Initial row index\r
 111    * @param j0\r
 112    *        Initial column index\r
 113    * @param nrows\r
 114    *        Number of rows\r
 115    * @param ncols\r
 116    *        Number of columns\r
 117    * @return submatrix\r
 118    * \r
 119    */\r
 120 \r
 121   public Matrix getSubmatrix(int i0, int j0, int nrows, int ncols) {\r
 122     Matrix x = new Matrix(null, nrows, ncols);\r
 123     double[][] xa = x.a;\r
 124     for (int i = nrows; --i >= 0;)\r
 125       for (int j = ncols; --j >= 0;)\r
 126         xa[i][j] = a[i0 + i][j0 + j];\r
 127     return x;\r
 128   }\r
 129 \r
 130   /**\r
 131    * Get a submatrix for a give number of columns and selected row set.\r
 132    * \r
 133    * @param r\r
 134    *        Array of row indices.\r
 135    * @param n\r
 136    *        number of rows \r
 137    * @return submatrix\r
 138    */\r
 139 \r
 140   public Matrix getMatrixSelected(int[] r, int n) {\r
 141     Matrix x = new Matrix(null, r.length, n);\r
 142     double[][] xa = x.a;\r
 143     for (int i = r.length; --i >= 0;) {\r
 144       double[] b = a[r[i]];\r
 145       for (int j = n; --j >= 0;)\r
 146         xa[i][j] = b[j];\r
 147     }\r
 148     return x;\r
 149   }\r
 150 \r
 151   /**\r
 152    * Matrix transpose.\r
 153    * \r
 154    * @return A'\r
 155    */\r
 156 \r
 157   public Matrix transpose() {\r
 158     Matrix x = new Matrix(null, n, m);\r
 159     double[][] c = x.a;\r
 160     for (int i = m; --i >= 0;)\r
 161       for (int j = n; --j >= 0;)\r
 162         c[j][i] = a[i][j];\r
 163     return x;\r
 164   }\r
 165 \r
 166   /**\r
 167    * add two matrices\r
 168    * @param b\r
 169    * @return new Matrix this + b\r
 170    */\r
 171   public Matrix add(Matrix b) {\r
 172     return scaleAdd(b, 1);\r
 173   }\r
 174 \r
 175   /**\r
 176    * subtract two matrices\r
 177    * @param b\r
 178    * @return new Matrix this - b\r
 179    */\r
 180   public Matrix sub(Matrix b) {\r
 181     return scaleAdd(b, -1);\r
 182   }\r
 183   \r
 184   /**\r
 185    * X = A + B*scale\r
 186    * @param b \r
 187    * @param scale \r
 188    * @return X\r
 189    * \r
 190    */\r
 191   public Matrix scaleAdd(Matrix b, double scale) {\r
 192     Matrix x = new Matrix(null, m, n);\r
 193     double[][] xa = x.a;\r
 194     double[][] ba = b.a;\r
 195     for (int i = m; --i >= 0;)\r
 196       for (int j = n; --j >= 0;)\r
 197         xa[i][j] = ba[i][j] * scale + a[i][j];\r
 198     return x;\r
 199   }\r
 200 \r
 201   /**\r
 202    * Linear algebraic matrix multiplication, A * B\r
 203    * \r
 204    * @param b\r
 205    *        another matrix\r
 206    * @return Matrix product, A * B or null for wrong dimension\r
 207    */\r
 208 \r
 209   public Matrix mul(Matrix b) {\r
 210     if (b.m != n)\r
 211       return null;\r
 212     Matrix x = new Matrix(null, m, b.n);\r
 213     double[][] xa = x.a;\r
 214     double[][] ba = b.a;\r
 215     for (int j = b.n; --j >= 0;)\r
 216       for (int i = m; --i >= 0;) {\r
 217         double[] arowi = a[i];\r
 218         double s = 0;\r
 219         for (int k = n; --k >= 0;)\r
 220           s += arowi[k] * ba[k][j];\r
 221         xa[i][j] = s;\r
 222       }\r
 223     return x;\r
 224   }\r
 225 \r
 226   /**\r
 227    * Matrix inverse or pseudoinverse\r
 228    * \r
 229    * @return inverse (m == n) or pseudoinverse (m != n)\r
 230    */\r
 231 \r
 232   public Matrix inverse() {\r
 233     return new LUDecomp(m, n).solve(identity(m, m), n);\r
 234   }\r
 235 \r
 236   /**\r
 237    * Matrix trace.\r
 238    * \r
 239    * @return sum of the diagonal elements.\r
 240    */\r
 241 \r
 242   public double trace() {\r
 243     double t = 0;\r
 244     for (int i = Math.min(m, n); --i >= 0;)\r
 245       t += a[i][i];\r
 246     return t;\r
 247   }\r
 248 \r
 249   /**\r
 250    * Generate identity matrix\r
 251    * \r
 252    * @param m\r
 253    *        Number of rows.\r
 254    * @param n\r
 255    *        Number of columns.\r
 256    * @return An m-by-n matrix with ones on the diagonal and zeros elsewhere.\r
 257    */\r
 258 \r
 259   public static Matrix identity(int m, int n) {\r
 260     Matrix x = new Matrix(null, m, n);\r
 261     double[][] xa = x.a;\r
 262     for (int i = Math.min(m, n); --i >= 0;)\r
 263       xa[i][i] = 1;\r
 264     return x;\r
 265   }\r
 266 \r
 267   /**\r
 268    * similarly to M3/M4 standard rotation/translation matrix\r
 269    * we set a rotationTranslation matrix to be:\r
 270    * \r
 271    * [   nxn rot    nx1 trans\r
 272    * \r
 273    *     1xn  0     1x1 1      ]\r
 274    * \r
 275    * \r
 276    * @return rotation matrix\r
 277    */\r
 278   public Matrix getRotation() {\r
 279     return getSubmatrix(0, 0, m - 1, n - 1);\r
 280   }\r
 281 \r
 282   public Matrix getTranslation() {\r
 283     return getSubmatrix(0, n - 1, m - 1, 1);\r
 284   }\r
 285 \r
 286   public static Matrix newT(T3 r, boolean asColumn) {\r
 287     return (asColumn ? new Matrix(new double[][] { new double[] { r.x },\r
 288         new double[] { r.y }, new double[] { r.z } }, 3, 1) : new Matrix(\r
 289         new double[][] { new double[] { r.x, r.y, r.z } }, 1, 3));\r
 290   }\r
 291 \r
 292   @Override\r
 293   public String toString() {\r
 294     String s = "[\n";\r
 295     for (int i = 0; i < m; i++) {\r
 296       s += "  [";\r
 297       for (int j = 0; j < n; j++)\r
 298         s += " " + a[i][j];\r
 299       s += "]\n";\r
 300     }\r
 301     s += "]";\r
 302     return s;\r
 303   }\r
 304 \r
 305   /**\r
 306    * \r
 307    * Edited down by Bob Hanson for minimum needed by Jmol -- just constructor\r
 308    * and solve\r
 309    * \r
 310    * LU Decomposition.\r
 311    * <P>\r
 312    * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit\r
 313    * lower triangular matrix L, an n-by-n upper triangular matrix U, and a\r
 314    * permutation vector piv of length m so that A(piv,:) = L*U. If m < n, then L\r
 315    * is m-by-m and U is m-by-n.\r
 316    * <P>\r
 317    * The LU decompostion with pivoting always exists, even if the matrix is\r
 318    * singular, so the constructor will never fail. The primary use of the LU\r
 319    * decomposition is in the solution of square systems of simultaneous linear\r
 320    * equations. This will fail if isNonsingular() returns false.\r
 321    */\r
 322 \r
 323   private class LUDecomp {\r
 324 \r
 325     /* ------------------------\r
 326        Class variables\r
 327      * ------------------------ */\r
 328 \r
 329     /**\r
 330      * Array for internal storage of decomposition.\r
 331      * \r
 332      */\r
 333     private double[][] LU;\r
 334 \r
 335     /**\r
 336      * Internal storage of pivot vector.\r
 337      * \r
 338      */\r
 339     private int[] piv;\r
 340 \r
 341     private int pivsign;\r
 342 \r
 343     /* ------------------------\r
 344        Constructor\r
 345      * ------------------------ */\r
 346 \r
 347     /**\r
 348      * LU Decomposition Structure to access L, U and piv.\r
 349      * @param m \r
 350      * @param n \r
 351      * \r
 352      */\r
 353 \r
 354     protected LUDecomp(int m, int n) {\r
 355       \r
 356       // Use a "left-looking", dot-product, Crout/Doolittle algorithm.\r
 357 \r
 358       LU = getArrayCopy();\r
 359       piv = new int[m];\r
 360       for (int i = m; --i >= 0;)\r
 361         piv[i] = i;\r
 362       pivsign = 1;\r
 363       double[] LUrowi;\r
 364       double[] LUcolj = new double[m];\r
 365 \r
 366       // Outer loop.\r
 367 \r
 368       for (int j = 0; j < n; j++) {\r
 369 \r
 370         // Make a copy of the j-th column to localize references.\r
 371 \r
 372         for (int i = m; --i >= 0;)\r
 373           LUcolj[i] = LU[i][j];\r
 374 \r
 375         // Apply previous transformations.\r
 376 \r
 377         for (int i = m; --i >= 0;) {\r
 378           LUrowi = LU[i];\r
 379 \r
 380           // Most of the time is spent in the following dot product.\r
 381 \r
 382           int kmax = Math.min(i, j);\r
 383           double s = 0.0;\r
 384           for (int k = kmax; --k >= 0;)\r
 385             s += LUrowi[k] * LUcolj[k];\r
 386 \r
 387           LUrowi[j] = LUcolj[i] -= s;\r
 388         }\r
 389 \r
 390         // Find pivot and exchange if necessary.\r
 391 \r
 392         int p = j;\r
 393         for (int i = m; --i > j;)\r
 394           if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p]))\r
 395             p = i;\r
 396         if (p != j) {\r
 397           for (int k = n; --k >= 0;) {\r
 398             double t = LU[p][k];\r
 399             LU[p][k] = LU[j][k];\r
 400             LU[j][k] = t;\r
 401           }\r
 402           int k = piv[p];\r
 403           piv[p] = piv[j];\r
 404           piv[j] = k;\r
 405           pivsign = -pivsign;\r
 406         }\r
 407 \r
 408         // Compute multipliers.\r
 409 \r
 410         if (j < m & LU[j][j] != 0.0)\r
 411           for (int i = m; --i > j;)\r
 412             LU[i][j] /= LU[j][j];\r
 413       }\r
 414     }\r
 415 \r
 416     /* ------------------------\r
 417        default Methods\r
 418      * ------------------------ */\r
 419 \r
 420     /**\r
 421      * Solve A*X = B\r
 422      * \r
 423      * @param b\r
 424      *        A Matrix with as many rows as A and any number of columns.\r
 425      * @param n \r
 426      * @return X so that L*U*X = B(piv,:) or null for wrong size or singular matrix\r
 427      */\r
 428 \r
 429     protected Matrix solve(Matrix b, int n) {\r
 430       for (int j = 0; j < n; j++)\r
 431         if (LU[j][j] == 0)\r
 432           return null; // matrix is singular\r
 433 \r
 434       // Copy right hand side with pivoting\r
 435       int nx = b.n;\r
 436       Matrix x = b.getMatrixSelected(piv, nx);\r
 437       double[][] a = x.a;\r
 438 \r
 439       // Solve L*Y = B(piv,:)\r
 440       for (int k = 0; k < n; k++)\r
 441         for (int i = k + 1; i < n; i++)\r
 442           for (int j = 0; j < nx; j++)\r
 443             a[i][j] -= a[k][j] * LU[i][k];\r
 444 \r
 445       // Solve U*X = Y;\r
 446       for (int k = n; --k >= 0;) {\r
 447         for (int j = nx; --j >= 0;)\r
 448           a[k][j] /= LU[k][k];\r
 449         for (int i = k; --i >= 0;)\r
 450           for (int j = nx; --j >= 0;)\r
 451             a[i][j] -= a[k][j] * LU[i][k];\r
 452       }\r
 453       return x;\r
 454     }\r
 455   }\r
 456 \r
 457 }\r