srcjar/javajs/util/Matrix.java

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