View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.numerics.math;
39  
40  /**
41   * The MatrixMath class is a simple matrix math library used mainly by the X-ray package.
42   *
43   * <p>All methods are thread-safe and static.
44   *
45   * @author Michael J. Schnieders
46   * @since 1.0
47   */
48  public final class MatrixMath {
49  
50    /**
51     * Returns the determinant for a 3x3 matrix.
52     *
53     * @param m input matrix.
54     * @return Returns the determinant.
55     */
56    public static double determinant3(double[][] m) {
57      return m[0][0] * m[1][1] * m[2][2]
58          - m[0][0] * m[1][2] * m[2][1]
59          + m[0][1] * m[1][2] * m[2][0]
60          - m[0][1] * m[1][0] * m[2][2]
61          + m[0][2] * m[1][0] * m[2][1]
62          - m[0][2] * m[1][1] * m[2][0];
63    }
64  
65    /**
66     * determinant3
67     *
68     * @param m an array of double.
69     * @return Returns the determinant.
70     */
71    public static double determinant3(double[] m) {
72      return m[0] * m[1] * m[2]
73          - m[0] * m[5] * m[5]
74          + m[3] * m[5] * m[4]
75          - m[3] * m[3] * m[2]
76          + m[4] * m[3] * m[5]
77          - m[4] * m[1] * m[4];
78    }
79  
80    /**
81     * inverse of a 3x3 matrix.
82     *
83     * @param m input matrix.
84     * @return Returns the matrix inverse.
85     */
86    public static double[][] mat3Inverse(double[][] m) {
87      return mat3Inverse(m, new double[3][3]);
88    }
89  
90    /**
91     * mat3inverse
92     *
93     * @param m an array of double.
94     * @param res an array of double.
95     * @return Returns the matrix res.
96     */
97    public static double[][] mat3Inverse(double[][] m, double[][] res) {
98      double det = determinant3(m);
99      res[0][0] = (m[1][1] * m[2][2] - m[1][2] * m[2][1]) / det;
100     res[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) / det;
101     res[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) / det;
102     res[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) / det;
103     res[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) / det;
104     res[1][2] = (m[0][2] * m[1][0] - m[0][0] * m[1][2]) / det;
105     res[2][0] = (m[1][0] * m[2][1] - m[1][1] * m[2][0]) / det;
106     res[2][1] = (m[0][1] * m[2][0] - m[0][0] * m[2][1]) / det;
107     res[2][2] = (m[0][0] * m[1][1] - m[0][1] * m[1][0]) / det;
108     return res;
109   }
110 
111   /**
112    * Matrix times a matrix.
113    *
114    * @param m1 First input matrix.
115    * @param m2 Second input matrix.
116    * @return Returns the matrix product.
117    */
118   public static double[][] mat3Mat3(double[][] m1, double[][] m2) {
119     return mat3Mat3(m1, m2, new double[3][3]);
120   }
121 
122   /**
123    * mat3mat3
124    *
125    * @param m1 an array of double.
126    * @param m2 an array of double.
127    * @param res an array of double.
128    * @return Returns the matrix res.
129    */
130   public static double[][] mat3Mat3(double[][] m1, double[][] m2, double[][] res) {
131     res[0][0] = m1[0][0] * m2[0][0] + m1[0][1] * m2[1][0] + m1[0][2] * m2[2][0];
132     res[0][1] = m1[0][0] * m2[0][1] + m1[0][1] * m2[1][1] + m1[0][2] * m2[2][1];
133     res[0][2] = m1[0][0] * m2[0][2] + m1[0][1] * m2[1][2] + m1[0][2] * m2[2][2];
134     res[1][0] = m1[1][0] * m2[0][0] + m1[1][1] * m2[1][0] + m1[1][2] * m2[2][0];
135     res[1][1] = m1[1][0] * m2[0][1] + m1[1][1] * m2[1][1] + m1[1][2] * m2[2][1];
136     res[1][2] = m1[1][0] * m2[0][2] + m1[1][1] * m2[1][2] + m1[1][2] * m2[2][2];
137     res[2][0] = m1[2][0] * m2[0][0] + m1[2][1] * m2[1][0] + m1[2][2] * m2[2][0];
138     res[2][1] = m1[2][0] * m2[0][1] + m1[2][1] * m2[1][1] + m1[2][2] * m2[2][1];
139     res[2][2] = m1[2][0] * m2[0][2] + m1[2][1] * m2[1][2] + m1[2][2] * m2[2][2];
140     return res;
141   }
142 
143   /**
144    * Matrix times a matrix (both 4x4).
145    *
146    * @param m1 First input matrix.
147    * @param m2 Second input matrix.
148    * @return Returns the matrix product.
149    */
150   public static double[][] mat4Mat4(double[][] m1, double[][] m2) {
151     return mat4Mat4(m1, m2, new double[4][4]);
152   }
153 
154   /**
155    * Multiply two 4x4 matrices.
156    *
157    * @param m1 an array of double (first matrix).
158    * @param m2 an array of double (second matrix).
159    * @param res Resultant matrix.
160    * @return Returns the matrix res.
161    */
162   public static double[][] mat4Mat4(double[][] m1, double[][] m2, double[][] res) {
163     assert (m1.length == 4);
164     assert (m1[0].length == 4);
165     assert (m2.length == 4);
166     assert (m2[0].length == 4);
167     assert (res.length == 4);
168     assert (res[0].length == 4);
169     res[0][0] =
170         m1[0][0] * m2[0][0] + m1[0][1] * m2[1][0] + m1[0][2] * m2[2][0] + m1[0][3] * m2[3][0];
171     res[0][1] =
172         m1[0][0] * m2[0][1] + m1[0][1] * m2[1][1] + m1[0][2] * m2[2][1] + m1[0][3] * m2[3][1];
173     res[0][2] =
174         m1[0][0] * m2[0][2] + m1[0][1] * m2[1][2] + m1[0][2] * m2[2][2] + m1[0][3] * m2[3][2];
175     res[0][3] =
176         m1[0][0] * m2[0][3] + m1[0][1] * m2[1][3] + m1[0][2] * m2[2][3] + m1[0][3] * m2[3][3];
177     res[1][0] =
178         m1[1][0] * m2[0][0] + m1[1][1] * m2[1][0] + m1[1][2] * m2[2][0] + m1[1][3] * m2[3][0];
179     res[1][1] =
180         m1[1][0] * m2[0][1] + m1[1][1] * m2[1][1] + m1[1][2] * m2[2][1] + m1[1][3] * m2[3][1];
181     res[1][2] =
182         m1[1][0] * m2[0][2] + m1[1][1] * m2[1][2] + m1[1][2] * m2[2][2] + m1[1][3] * m2[3][2];
183     res[1][3] =
184         m1[1][0] * m2[0][3] + m1[1][1] * m2[1][3] + m1[1][2] * m2[2][3] + m1[1][3] * m2[3][3];
185     res[2][0] =
186         m1[2][0] * m2[0][0] + m1[2][1] * m2[1][0] + m1[2][2] * m2[2][0] + m1[2][3] * m2[3][0];
187     res[2][1] =
188         m1[2][0] * m2[0][1] + m1[2][1] * m2[1][1] + m1[2][2] * m2[2][1] + m1[2][3] * m2[3][1];
189     res[2][2] =
190         m1[2][0] * m2[0][2] + m1[2][1] * m2[1][2] + m1[2][2] * m2[2][2] + m1[2][3] * m2[3][2];
191     res[2][3] =
192         m1[2][0] * m2[0][3] + m1[2][1] * m2[1][3] + m1[2][2] * m2[2][3] + m1[2][3] * m2[3][3];
193     res[3][0] =
194         m1[3][0] * m2[0][0] + m1[3][1] * m2[1][0] + m1[3][2] * m2[2][0] + m1[3][3] * m2[3][0];
195     res[3][1] =
196         m1[3][0] * m2[0][1] + m1[3][1] * m2[1][1] + m1[3][2] * m2[2][1] + m1[3][3] * m2[3][1];
197     res[3][2] =
198         m1[3][0] * m2[0][2] + m1[3][1] * m2[1][2] + m1[3][2] * m2[2][2] + m1[3][3] * m2[3][2];
199     res[3][3] =
200         m1[3][0] * m2[0][3] + m1[3][1] * m2[1][3] + m1[3][2] * m2[2][3] + m1[3][3] * m2[3][3];
201     return res;
202   }
203 
204   /**
205    * matrix times a vector representation of a symmetric 3x3 matrix
206    *
207    * @param m input matrix
208    * @param v input vector of the form 11, 22, 33, 12, 13, 23
209    * @return matrix product
210    */
211   public static double[][] mat3SymVec6(double[][] m, double[] v) {
212     return mat3SymVec6(m, v, new double[3][3]);
213   }
214 
215   /**
216    * mat3SymVec6
217    *
218    * @param m an array of double.
219    * @param v an array of double.
220    * @param res an array of double.
221    * @return Returns the matrix res.
222    */
223   public static double[][] mat3SymVec6(double[][] m, double[] v, double[][] res) {
224     res[0][0] = m[0][0] * v[0] + m[0][1] * v[3] + m[0][2] * v[4];
225     res[0][1] = m[0][0] * v[3] + m[0][1] * v[1] + m[0][2] * v[5];
226     res[0][2] = m[0][0] * v[4] + m[0][1] * v[5] + m[0][2] * v[2];
227     res[1][0] = m[1][0] * v[0] + m[1][1] * v[3] + m[1][2] * v[4];
228     res[1][1] = m[1][0] * v[3] + m[1][1] * v[1] + m[1][2] * v[5];
229     res[1][2] = m[1][0] * v[4] + m[1][1] * v[5] + m[1][2] * v[2];
230     res[2][0] = m[2][0] * v[0] + m[2][1] * v[3] + m[2][2] * v[4];
231     res[2][1] = m[2][0] * v[3] + m[2][1] * v[1] + m[2][2] * v[5];
232     res[2][2] = m[2][0] * v[4] + m[2][1] * v[5] + m[2][2] * v[2];
233     return res;
234   }
235 
236   /**
237    * matrix times a vector
238    *
239    * @param m input matrix.
240    * @param v input vector.
241    * @return Returns the vector product.
242    */
243   public static double[] mat3Vec3(double[] v, double[][] m) {
244     return mat3Vec3(v, m, new double[3]);
245   }
246 
247   /**
248    * mat3vec3
249    *
250    * @param v an array of double.
251    * @param m an array of double.
252    * @param res an array of double.
253    * @return Returns the matrix res.
254    */
255   public static double[] mat3Vec3(double[] v, double[][] m, double[] res) {
256     res[0] = m[0][0] * v[0] + m[0][1] * v[1] + m[0][2] * v[2];
257     res[1] = m[1][0] * v[0] + m[1][1] * v[1] + m[1][2] * v[2];
258     res[2] = m[2][0] * v[0] + m[2][1] * v[1] + m[2][2] * v[2];
259     return res;
260   }
261 
262   /**
263    * scalar times a matrix times a matrix.
264    *
265    * @param scalar Input scalar.
266    * @param m1 First input matrix.
267    * @param m2 Second input matrix.
268    * @return Returns the matrix product.
269    */
270   public static double[][] scalarMat3Mat3(double scalar, double[][] m1, double[][] m2) {
271     return scalarMat3Mat3(scalar, m1, m2, new double[3][3]);
272   }
273 
274   /**
275    * scalarMat3mat3
276    *
277    * @param scalar a double.
278    * @param m1 an array of double.
279    * @param m2 an array of double.
280    * @param res an array of double.
281    * @return Returns the matrix res.
282    */
283   public static double[][] scalarMat3Mat3(
284       double scalar, double[][] m1, double[][] m2, double[][] res) {
285     res[0][0] =
286         (scalar * m1[0][0]) * m2[0][0]
287             + (scalar * m1[0][1]) * m2[1][0]
288             + (scalar * m1[0][2]) * m2[2][0];
289     res[0][1] =
290         (scalar * m1[0][0]) * m2[0][1]
291             + (scalar * m1[0][1]) * m2[1][1]
292             + (scalar * m1[0][2]) * m2[2][1];
293     res[0][2] =
294         (scalar * m1[0][0]) * m2[0][2]
295             + (scalar * m1[0][1]) * m2[1][2]
296             + (scalar * m1[0][2]) * m2[2][2];
297     res[1][0] =
298         (scalar * m1[1][0]) * m2[0][0]
299             + (scalar * m1[1][1]) * m2[1][0]
300             + (scalar * m1[1][2]) * m2[2][0];
301     res[1][1] =
302         (scalar * m1[1][0]) * m2[0][1]
303             + (scalar * m1[1][1]) * m2[1][1]
304             + (scalar * m1[1][2]) * m2[2][1];
305     res[1][2] =
306         (scalar * m1[1][0]) * m2[0][2]
307             + (scalar * m1[1][1]) * m2[1][2]
308             + (scalar * m1[1][2]) * m2[2][2];
309     res[2][0] =
310         (scalar * m1[2][0]) * m2[0][0]
311             + (scalar * m1[2][1]) * m2[1][0]
312             + (scalar * m1[2][2]) * m2[2][0];
313     res[2][1] =
314         (scalar * m1[2][0]) * m2[0][1]
315             + (scalar * m1[2][1]) * m2[1][1]
316             + (scalar * m1[2][2]) * m2[2][1];
317     res[2][2] =
318         (scalar * m1[2][0]) * m2[0][2]
319             + (scalar * m1[2][1]) * m2[1][2]
320             + (scalar * m1[2][2]) * m2[2][2];
321     return res;
322   }
323 
324   /**
325    * vector representation of a symmetric 3x3 matrix times a matrix
326    *
327    * @param v input vector of the form 11, 22, 33, 12, 13, 23.
328    * @param m input matrix.
329    * @return Returns matrix product.
330    */
331   public static double[][] symVec6Mat3(double[] v, double[][] m) {
332     double[][] res = new double[3][3];
333     symVec6Mat3(v, m, res);
334     return res;
335   }
336 
337   /**
338    * symVec6mat3
339    *
340    * @param v an array of double.
341    * @param m an array of double.
342    * @param res an array of double.
343    * @return Returns the matrix res.
344    */
345   public static double[][] symVec6Mat3(double[] v, double[][] m, double[][] res) {
346     res[0][0] = v[0] * m[0][0] + v[3] * m[1][0] + v[4] * m[2][0];
347     res[0][1] = v[0] * m[0][1] + v[3] * m[1][1] + v[4] * m[2][1];
348     res[0][2] = v[0] * m[0][2] + v[3] * m[1][2] + v[4] * m[2][2];
349     res[1][0] = v[3] * m[0][0] + v[1] * m[1][0] + v[5] * m[2][0];
350     res[1][1] = v[3] * m[0][1] + v[1] * m[1][1] + v[5] * m[2][1];
351     res[1][2] = v[3] * m[0][2] + v[1] * m[1][2] + v[5] * m[2][2];
352     res[2][0] = v[4] * m[0][0] + v[5] * m[1][0] + v[2] * m[2][0];
353     res[2][1] = v[4] * m[0][1] + v[5] * m[1][1] + v[2] * m[2][1];
354     res[2][2] = v[4] * m[0][2] + v[5] * m[1][2] + v[2] * m[2][2];
355     return res;
356   }
357 
358   /**
359    * transpose3
360    *
361    * @param m Matrix m.
362    * @return Returns the transposed matrix in a new array.
363    */
364   public static double[][] transpose3(double[][] m) {
365     return transpose3(m, new double[3][3]);
366   }
367 
368   /**
369    * transpose3
370    *
371    * @param m Matrix m.
372    * @param t Matrix t.
373    * @return Returns the matrix t.
374    */
375   public static double[][] transpose3(double[][] m, double[][] t) {
376     t[0][0] = m[0][0];
377     t[0][1] = m[1][0];
378     t[0][2] = m[2][0];
379     t[1][0] = m[0][1];
380     t[1][1] = m[1][1];
381     t[1][2] = m[2][1];
382     t[2][0] = m[0][2];
383     t[2][1] = m[1][2];
384     t[2][2] = m[2][2];
385     return t;
386   }
387 
388   /**
389    * Vector times a matrix.
390    *
391    * @param v input vector.
392    * @param m input matrix.
393    * @return Returns the vector product.
394    */
395   public static double[] vec3Mat3(double[] v, double[][] m) {
396     return vec3Mat3(v, m, new double[3]);
397   }
398 
399   /**
400    * vec3mat3
401    *
402    * @param v an array of double.
403    * @param m an array of double.
404    * @param res an array of double.
405    * @return Returns the array res.
406    */
407   public static double[] vec3Mat3(double[] v, double[][] m, double[] res) {
408     res[0] = v[0] * m[0][0] + v[1] * m[1][0] + v[2] * m[2][0];
409     res[1] = v[0] * m[0][1] + v[1] * m[1][1] + v[2] * m[2][1];
410     res[2] = v[0] * m[0][2] + v[1] * m[1][2] + v[2] * m[2][2];
411     return res;
412   }
413 
414   /**
415    * Multiply coordinates by the transpose of a matrix.
416    *
417    * @param in input coordinates.
418    * @param out output coordinates.
419    * @param matrix multiply by the transpose of this matrix.
420    */
421   public static void applyMatrixTranspose(double[] in, double[] out, double[][] matrix) {
422     double xc = in[0];
423     double yc = in[1];
424     double zc = in[2];
425     out[0] = xc * matrix[0][0] + yc * matrix[1][0] + zc * matrix[2][0];
426     out[1] = xc * matrix[0][1] + yc * matrix[1][1] + zc * matrix[2][1];
427     out[2] = xc * matrix[0][2] + yc * matrix[1][2] + zc * matrix[2][2];
428   }
429 }