hoomd_linear_algebra/
matrix.rs

1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4use serde::{Deserialize, Serialize};
5use serde_with::serde_as;
6use std::fmt;
7
8use crate::{Full, GeneralMatrix, Invertible, MatMul, QuadraticForm, SquareMatrix};
9
10/// ``std::ops`` implementations for [`Matrix`]
11mod ops;
12
13pub use crate::diagonal::DiagonalMatrix;
14
15/// A matrix with N rows and M columns, allocated on the stack.
16///
17/// # Example
18/// ```
19/// use hoomd_linear_algebra::matrix::Matrix;
20///
21/// let a = Matrix {
22///     rows: [[-1.0, 4.0, -6.0], [2.0, -3.0, 1.0]],
23/// };
24/// ```
25#[serde_as]
26#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
27pub struct Matrix<const N: usize, const M: usize> {
28    /// The elements of the matrix
29    #[serde_as(as = "[[_; M]; N]")]
30    pub rows: [[f64; M]; N],
31}
32/// A 2x2 matrix, allocated on the stack.
33pub type Matrix22 = Matrix<2, 2>;
34/// A 3x3 matrix, allocated on the stack.
35pub type Matrix33 = Matrix<3, 3>;
36/// A 4x4 matrix, allocated on the stack.
37pub type Matrix44 = Matrix<4, 4>;
38
39/// cos(π/8), used in svd 3x3
40const C_STAR: f64 = 0.923_879_532_511_286_8;
41/// sin(π/8), used in svd 3x3
42const S_STAR: f64 = 0.382_683_432_365_089_8;
43
44impl<const N: usize, const M: usize> GeneralMatrix for Matrix<N, M> {
45    /// Fill a matrix with zeros.
46    ///
47    /// # Examples
48    /// ```
49    /// use hoomd_linear_algebra::{GeneralMatrix, matrix::Matrix22};
50    ///
51    /// let m = Matrix22::zeros();
52    /// assert_eq!(m.rows, [[0.0, 0.0], [0.0, 0.0]]);
53    /// ```
54    #[inline]
55    fn zeros() -> Self {
56        Self {
57            rows: std::array::from_fn(|_| std::array::from_fn(|_| 0.0)),
58        }
59    }
60    #[inline]
61    fn shape(&self) -> (usize, usize) {
62        (self.n_rows(), self.n_columns())
63    }
64}
65
66impl<const N: usize, const M: usize> Default for Matrix<N, M> {
67    /// The default matrix is filled with zeros.
68    ///
69    /// # Examples
70    /// ```
71    /// use hoomd_linear_algebra::matrix::Matrix22;
72    ///
73    /// let m = Matrix22::default();
74    /// assert_eq!(m.rows, [[0.0, 0.0], [0.0, 0.0]]);
75    /// ```
76    #[inline]
77    fn default() -> Self {
78        Self::zeros()
79    }
80}
81
82impl<const N: usize, const M: usize> Full for Matrix<N, M> {
83    /// Construct a matrix with the same value in every element.
84    ///
85    /// # Examples
86    /// ```
87    /// use hoomd_linear_algebra::{Full, matrix::Matrix22};
88    ///
89    /// let m = Matrix22::full(5.0);
90    /// assert_eq!(m.rows, [[5.0, 5.0], [5.0, 5.0]]);
91    /// ```
92    #[inline]
93    fn full(value: f64) -> Self {
94        Self {
95            rows: std::array::from_fn(|_| std::array::from_fn(|_| value)),
96        }
97    }
98}
99
100impl<const N: usize> SquareMatrix for Matrix<N, N> {
101    /// # Examples
102    /// ```
103    /// use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};
104    ///
105    /// let m = Matrix22::identity();
106    /// assert_eq!(m.rows, [[1.0, 0.0], [0.0, 1.0]]);
107    /// ```
108    #[inline]
109    fn identity() -> Self {
110        Self {
111            rows: std::array::from_fn(|i| std::array::from_fn(|j| if i == j { 1.0 } else { 0.0 })),
112        }
113    }
114}
115
116impl<const N: usize, const M: usize, const K: usize> MatMul<Matrix<M, K>> for Matrix<N, M> {
117    type Output = Matrix<N, K>;
118    /// Matrix-matrix multiplication.
119    ///
120    /// # Examples
121    /// ```
122    /// use hoomd_linear_algebra::{
123    ///     MatMul,
124    ///     matrix::{Matrix, Matrix22},
125    /// };
126    ///
127    /// let a = Matrix22 {
128    ///     rows: [[1.0, 2.0], [3.0, 4.0]],
129    /// };
130    /// let b = Matrix22 {
131    ///     rows: [[5.0, 6.0], [7.0, 8.0]],
132    /// };
133    /// let c = a.matmul(&b);
134    /// assert_eq!(c.rows, [[19.0, 22.0], [43.0, 50.0]]);
135    /// ```
136    #[inline]
137    fn matmul(&self, rhs: &Matrix<M, K>) -> Self::Output {
138        let mut result = Self::Output::zeros();
139        for n in 0..N {
140            for k in 0..K {
141                for m in 0..M {
142                    result.rows[n][k] += self.rows[n][m] * rhs.rows[m][k];
143                }
144            }
145        }
146
147        result
148    }
149}
150
151impl<const N: usize, const M: usize> MatMul<DiagonalMatrix<M>> for Matrix<N, M> {
152    type Output = Matrix<N, M>;
153
154    /// Matrix-diagonal matrix multiplication.
155    ///
156    /// This is equivalent to scaling each column of a [`Matrix`] by the corresponding
157    /// element in a [`DiagonalMatrix`].
158    ///
159    /// # Examples
160    /// ```
161    /// use hoomd_linear_algebra::{
162    ///     Full, GeneralMatrix, MatMul,
163    ///     matrix::{DiagonalMatrix, Matrix22},
164    /// };
165    ///
166    /// let diag = DiagonalMatrix {
167    ///     elements: [3.0, 4.0],
168    /// };
169    /// let a = Matrix22::full(1.0).matmul(&diag);
170    /// assert_eq!(a[(0, 1)], 4.0);
171    /// assert_eq!(a[(1, 0)], 3.0);
172    /// ```
173    #[inline]
174    fn matmul(&self, rhs: &DiagonalMatrix<M>) -> Self::Output {
175        let mut result = Self::Output::zeros();
176        for (i, row) in result.rows.iter_mut().enumerate().take(N) {
177            for j in 0..M {
178                row[j] = self.rows[i][j] * rhs[j];
179            }
180        }
181        result
182    }
183}
184
185impl<const N: usize, const M: usize> Matrix<N, M> {
186    /// Interchange the rows and columns of matrix.
187    ///
188    /// ```math
189    /// \mathbf{A}^\top_{ji} = \mathbf{A}_{ij}
190    /// ```
191    ///
192    /// # Examples
193    /// ```
194    /// use hoomd_linear_algebra::matrix::Matrix;
195    ///
196    /// let m = Matrix {
197    ///     rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
198    /// };
199    /// let m_t = m.transpose();
200    /// assert_eq!(m_t.rows, [[1.0, 4.0], [2.0, 5.0], [3.0, 6.0]]);
201    /// ```
202    #[inline]
203    #[must_use]
204    pub fn transpose(&self) -> Matrix<M, N> {
205        Matrix {
206            rows: std::array::from_fn(|j| std::array::from_fn(|i| self[(i, j)])),
207        }
208    }
209
210    /// Apply a function to an [`Matrix`] by rows.
211    ///
212    /// # Examples
213    /// ```
214    /// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};
215    ///
216    /// let m = Matrix22::full(3.0);
217    /// assert_eq!(
218    ///     m.map_rows(|v| [v[0] + 2.0, v[1]]),
219    ///     Matrix22 {
220    ///         rows: [[5.0, 3.0], [5.0, 3.0]]
221    ///     }
222    /// );
223    /// ```
224    #[inline]
225    #[must_use]
226    pub fn map_rows<F>(self, f: F) -> Self
227    where
228        F: FnMut([f64; M]) -> [f64; M],
229    {
230        Self {
231            rows: self.rows.map(f),
232        }
233    }
234
235    /// Apply a function to an [`Matrix`] by columns.
236    ///
237    /// # Examples
238    /// ```
239    /// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};
240    ///
241    /// let m = Matrix22::full(3.0);
242    /// assert_eq!(
243    ///     m.map_columns(|v| [v[0] + 2.0, v[1]]),
244    ///     Matrix22 {
245    ///         rows: [[5.0, 5.0], [3.0, 3.0]]
246    ///     }
247    /// );
248    /// ```
249    #[inline]
250    #[must_use]
251    pub fn map_columns<F>(self, f: F) -> Self
252    where
253        F: FnMut([f64; N]) -> [f64; N],
254    {
255        self.clone().transpose().map_rows(f).transpose()
256    }
257    /// Apply a function to [`Matrix`] elements.
258    ///
259    /// # Examples
260    /// ```
261    /// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix33};
262    ///
263    /// let m = Matrix33::full(3.0);
264    /// assert_eq!(m.map_elements(|x| x + 2.0), m + Matrix33::full(2.0));
265    /// ```
266    #[inline]
267    #[must_use]
268    pub fn map_elements<F>(self, f: F) -> Self
269    where
270        F: Fn(f64) -> f64,
271    {
272        Self {
273            rows: self.rows.map(|v| v.map(&f)),
274        }
275    }
276
277    /// Iterate over every element in the [`Matrix`].
278    ///
279    /// The iterator yields all matrix elements in row major order.
280    ///
281    /// # Examples
282    /// ```
283    /// use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};
284    ///
285    /// let x = Matrix22 {
286    ///     rows: [[1.0, 2.0], [3.0, 4.0]],
287    /// };
288    ///
289    /// let mut iterator = x.iter_elements();
290    /// assert_eq!(iterator.next(), Some(1.0));
291    /// assert_eq!(iterator.next(), Some(2.0));
292    /// assert_eq!(iterator.next(), Some(3.0));
293    /// assert_eq!(iterator.next(), Some(4.0));
294    /// assert_eq!(iterator.next(), None);
295    /// ```
296    #[inline]
297    pub fn iter_elements(&self) -> impl Iterator<Item = f64> {
298        self.rows.iter().flat_map(|row| row.iter().copied())
299    }
300
301    /// Iterate over every element in the [`Matrix`].
302    ///
303    /// The iterator yields all matrix elements in row major order.
304    ///
305    /// # Examples
306    /// ```
307    /// use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};
308    ///
309    /// let mut x = Matrix22 {
310    ///     rows: [[1.0, 2.0], [3.0, 4.0]],
311    /// };
312    ///
313    /// let x_copy = x.clone();
314    /// let mut iterator = x.iter_elements_mut();
315    /// iterator.for_each(|x| *x *= 2.0);
316    /// assert_eq!(x, x_copy * 2.0);
317    /// ```
318    #[inline]
319    pub fn iter_elements_mut(&mut self) -> impl Iterator<Item = &mut f64> {
320        self.rows.iter_mut().flat_map(|row| row.iter_mut())
321    }
322
323    /// Iterate over matrix rows.
324    ///
325    /// # Examples
326    /// ```
327    /// use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};
328    ///
329    /// let x = Matrix22 {
330    ///     rows: [[1.0, 2.0], [3.0, 4.0]],
331    /// };
332    /// let mut iterator = x.iter_rows();
333    /// assert_eq!(iterator.next(), Some([1.0, 2.0]));
334    /// assert_eq!(iterator.next(), Some([3.0, 4.0]));
335    /// assert_eq!(iterator.next(), None);
336    /// ```
337    #[inline]
338    pub fn iter_rows(&self) -> impl Iterator<Item = [f64; M]> {
339        self.rows.iter().copied()
340    }
341
342    /// Folds every element into an accumulator by applying an operation.
343    ///
344    /// `fold_elements` takes two arguments: an initial value, and a closure
345    /// with two arguments: an ‘accumulator’, and an element. The closure
346    /// returns the value that the accumulator should have for the next iteration.
347    ///
348    /// The initial value is the value the accumulator will have on the first call.
349    /// After applying this closure to every element of the flattened iterator,
350    /// `fold_elements` returns the accumulator.
351    ///
352    /// # Examples
353    /// ```
354    /// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix22};
355    ///
356    /// let m = Matrix22::full(3.0);
357    /// // Sum the elements of a matrix
358    /// assert_eq!(m.fold_elements(0.0, |acc, x| acc + x), 3.0 * 4.0);
359    /// ```
360    #[inline]
361    pub fn fold_elements<B, F>(self, init: B, mut f: F) -> B
362    where
363        F: FnMut(B, f64) -> B,
364    {
365        let mut accum = init;
366        for x in self.iter_elements() {
367            accum = f(accum, x);
368        }
369        accum
370    }
371
372    /// Folds every row into an accumulator by applying an operation.
373    ///
374    /// `fold_rows` takes two arguments: an initial value, and a closure with two
375    /// arguments: an ‘accumulator’, and an element. The closure returns the
376    /// value that the accumulator should have for the next iteration.
377    ///
378    /// The initial value is the value the accumulator will have on the first
379    /// call. After applying this closure to every element of the flattened
380    /// iterator, `fold_rows` returns the accumulator.
381    ///
382    /// # Examples
383    /// ```
384    /// use hoomd_linear_algebra::{GeneralMatrix, matrix::Matrix22};
385    ///
386    /// let m = Matrix22 {
387    ///     rows: [[1.0, 2.0], [3.0, 4.0]],
388    /// };
389    /// // Average the columns of a matrix
390    /// let n_rows = m.n_rows() as f64;
391    /// assert_eq!(
392    ///     m.fold_rows([0.0; 2], |acc, x| [acc[0] + x[0], acc[1] + x[1]])
393    ///         .map(|x| x / n_rows),
394    ///     [(1.0 + 3.0) / 2.0, (2.0 + 4.0) / 2.0]
395    /// );
396    /// ```
397    #[inline]
398    pub fn fold_rows<B, F>(self, init: B, mut f: F) -> B
399    where
400        F: FnMut(B, [f64; M]) -> B,
401    {
402        let mut accum = init;
403        for x in self.iter_rows() {
404            accum = f(accum, x);
405        }
406        accum
407    }
408
409    /// Get the number of rows in the [`Matrix`].
410    ///
411    /// # Examples
412    /// ```
413    /// use hoomd_linear_algebra::matrix::Matrix;
414    ///
415    /// let m = Matrix {
416    ///     rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
417    /// };
418    /// assert_eq!(m.n_rows(), 2);
419    /// ```
420    #[must_use]
421    #[inline]
422    pub const fn n_rows(&self) -> usize {
423        N
424    }
425    /// Get the number of columns in the [`Matrix`].
426    ///
427    /// # Examples
428    /// ```
429    /// use hoomd_linear_algebra::matrix::Matrix;
430    ///
431    /// let m = Matrix {
432    ///     rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
433    /// };
434    /// assert_eq!(m.n_columns(), 3);
435    /// ```
436    #[must_use]
437    #[inline]
438    pub const fn n_columns(&self) -> usize {
439        M
440    }
441}
442impl<const N: usize> Matrix<N, N> {
443    /// Construct a square matrix with the given diagonal.
444    ///
445    /// `diagonal[i]` is the matrix element $` A_{ii} `$.
446    /// All off-diagonal elements are 0.
447    ///
448    /// # Example
449    ///
450    /// ```
451    /// use hoomd_linear_algebra::matrix::Matrix;
452    ///
453    /// let a = Matrix::with_diagonal([2.0, -3.0]);
454    ///
455    /// assert_eq!(a.rows, [[2.0, 0.0], [0.0, -3.0]]);
456    /// ```
457    #[inline]
458    #[must_use]
459    pub fn with_diagonal(diagonal: [f64; N]) -> Self {
460        DiagonalMatrix { elements: diagonal }.to_dense()
461    }
462
463    /// Compute the signed hypervolume of the hyperparallelepiped defined by a matrix.
464    ///
465    /// This implementation uses the Laplace expansion, which is optimal for small
466    /// matrices but will be extremely slow for large matrices due to its $`O(N!)`$
467    /// complexity.
468    ///
469    /// # Example
470    ///
471    /// ```
472    /// use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};
473    ///
474    /// let identity = Matrix22::identity();
475    /// assert_eq!(identity.determinant(), 1.0);
476    ///
477    /// let scaled = identity * 2.0;
478    /// assert_eq!(scaled.determinant(), 2.0 * 2.0);
479    /// ```
480    #[must_use]
481    #[inline]
482    pub fn determinant(&self) -> f64 {
483        // Compute the determinant of a 2x2 minor.
484        #[inline]
485        fn det2(a: f64, b: f64, c: f64, d: f64) -> f64 {
486            a * d - b * c
487        }
488        // Because math with const generics is not allowed in rust, we compute the indices
489        // of each submatrix and recur on those noncontiguous segments of the input.
490        #[inline]
491        fn det_recursive_noslice<const N: usize>(
492            matrix: &Matrix<N, N>,
493            row: usize,
494            col_indices: [usize; N],
495            minor_size: usize,
496        ) -> f64 {
497            // If we recur any lower than 4x4 minors, performance drops dramatically
498            if minor_size == 4 {
499                let r = matrix.rows;
500                let c = col_indices;
501
502                // Map recursive indices to direct matrix indices
503                let (i0, i1, i2, i3) = (row, row + 1, row + 2, row + 3);
504                let [j0, j1, j2, j3] = c[..4] else {
505                    unreachable!() // N >= 4 if we reach this point
506                };
507
508                let m0 = det2(r[i2][j2], r[i2][j3], r[i3][j2], r[i3][j3]);
509                let m1 = det2(r[i2][j1], r[i2][j3], r[i3][j1], r[i3][j3]);
510                let m2 = det2(r[i2][j1], r[i2][j2], r[i3][j1], r[i3][j2]);
511                let m3 = det2(r[i2][j0], r[i2][j3], r[i3][j0], r[i3][j3]);
512                let m4 = det2(r[i2][j0], r[i2][j2], r[i3][j0], r[i3][j2]);
513                let m5 = det2(r[i2][j0], r[i2][j1], r[i3][j0], r[i3][j1]);
514
515                return r[i0][j0] * (r[i1][j1] * m0 - r[i1][j2] * m1 + r[i1][j3] * m2)
516                    - r[i0][j1] * (r[i1][j0] * m0 - r[i1][j2] * m3 + r[i1][j3] * m4)
517                    + r[i0][j2] * (r[i1][j0] * m1 - r[i1][j1] * m3 + r[i1][j3] * m5)
518                    - r[i0][j3] * (r[i1][j0] * m2 - r[i1][j1] * m4 + r[i1][j2] * m5);
519            }
520
521            (0..minor_size).fold(0.0, |acc, idx| {
522                let minor_size = minor_size - 1;
523                let mut minor_cols = [0; N];
524                for j in 0..minor_size {
525                    // Store the indices for the next recursion, skipping col idx
526                    minor_cols[j] = col_indices[j + usize::from(j >= idx)];
527                }
528
529                let sign = if idx % 2 == 0 { 1.0 } else { -1.0 };
530                acc + sign
531                    * matrix.rows[row][col_indices[idx]]
532                    * det_recursive_noslice(matrix, row + 1, minor_cols, minor_size)
533            })
534        }
535        // Early exit for small matrices to ensure we get the optimal code.
536        match N {
537            0 => return 0.0,
538            1 => return self[(0, 0)],
539            2 => return det2(self[(0, 0)], self[(1, 0)], self[(0, 1)], self[(1, 1)]),
540            3 => {
541                return self[(0, 0)] * det2(self[(1, 1)], self[(1, 2)], self[(2, 1)], self[(2, 2)])
542                    - self[(0, 1)] * det2(self[(1, 0)], self[(1, 2)], self[(2, 0)], self[(2, 2)])
543                    + self[(0, 2)] * det2(self[(1, 0)], self[(1, 1)], self[(2, 0)], self[(2, 1)]);
544            }
545            _ => (),
546        }
547
548        let col_indices = std::array::from_fn(|i| i);
549        det_recursive_noslice(self, 0, col_indices, N)
550    }
551
552    /// Compute the sum of diagonal elements of a square matrix.
553    ///
554    /// # Examples
555    ///
556    /// ```
557    /// use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};
558    ///
559    /// let identity = Matrix22::identity();
560    /// assert_eq!(identity.trace(), 2.0);
561    ///
562    /// let scaled = identity * 3.0;
563    /// assert_eq!(scaled.trace(), 3.0 + 3.0);
564    /// ```
565    #[must_use]
566    #[inline]
567    pub fn trace(&self) -> f64 {
568        std::array::from_fn::<_, N, _>(|i| self[(i, i)])
569            .iter()
570            .sum()
571    }
572
573    /// Compute a matrix to an integer power
574    ///
575    /// ```math
576    /// \mathbf{A}^n = \prod_{i=1}^n \mathbf{A}
577    /// ```
578    ///
579    /// # Examples
580    ///
581    /// ```
582    /// use hoomd_linear_algebra::{Full, GeneralMatrix, MatMul, matrix::Matrix22};
583    ///
584    /// let matrix = Matrix22::full(2.0);
585    ///
586    /// assert_eq!(matrix.powi(2), matrix.matmul(&matrix));
587    ///
588    /// assert_eq!(matrix.powi(2).powi(2), matrix.powi(4));
589    /// ```
590    #[must_use]
591    #[inline]
592    pub fn powi(&self, n: u32) -> Self {
593        (0..n).fold(Self::identity(), |acc, _| acc.matmul(self))
594    }
595
596    /// Extract the diagonal elements from a square matrix.
597    ///
598    /// This method returns a `DiagonalMatrix<N>` containing the diagonal elements
599    /// of the input matrix, where the element at position `(i, i)` is taken from
600    /// the input matrix. All off-diagonal elements are ignored.
601    ///
602    /// # Examples
603    /// ```
604    /// use hoomd_linear_algebra::matrix::Matrix33;
605    /// let a = Matrix33 {
606    ///     rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]],
607    /// };
608    /// let b = a.diagonal();
609    /// assert_eq!(b.elements, [1.0, 5.0, 9.0]);
610    /// ```
611    #[must_use]
612    #[inline]
613    pub fn diagonal(&self) -> DiagonalMatrix<N> {
614        DiagonalMatrix {
615            elements: std::array::from_fn(|i| self.rows[i][i]),
616        }
617    }
618}
619
620impl<const N: usize> QuadraticForm<N> for Matrix<N, N> {
621    #[inline]
622    fn compute_quadratic_form(&self, x: &[f64; N]) -> f64 {
623        let mut result = 0.0;
624
625        for i in 0..N {
626            for j in 0..N {
627                result += x[i] * self[(i, j)] * x[j];
628            }
629        }
630        result
631    }
632}
633
634impl Invertible for Matrix<2, 2> {
635    /// Compute the inverse of a matrix. Will be `None` if the matrix is not invertible.
636    ///
637    /// This implementation uses a closed form solution for the matrix inverse.
638    ///
639    /// # Examples
640    /// ```
641    /// use hoomd_linear_algebra::{Invertible, matrix::Matrix22};
642    ///
643    /// let m = Matrix22 {
644    ///     rows: [[1.0, 2.0], [3.0, 4.0]],
645    /// };
646    /// let m_inv = m.inverse().unwrap();
647    /// assert_eq!(m_inv.rows, [[-2.0, 1.0], [1.5, -0.5]]);
648    /// ```
649    #[inline]
650    fn inverse(&self) -> Option<Self> {
651        let det = self.determinant();
652        if det == 0.0 {
653            None
654        } else {
655            let inv_det = det.recip();
656            Some(Self {
657                rows: [
658                    [inv_det * self.rows[1][1], inv_det * -self.rows[0][1]],
659                    [inv_det * -self.rows[1][0], inv_det * self.rows[0][0]],
660                ],
661            })
662        }
663    }
664}
665
666impl Invertible for Matrix<3, 3> {
667    /// Compute the inverse of a matrix. Will be `None` if the matrix is not invertible.
668    ///
669    /// This implementation uses a closed form solution for the matrix inverse based on
670    /// the cross product of rows.
671    ///
672    /// # Example
673    /// ```
674    /// use hoomd_linear_algebra::{Invertible, SquareMatrix, matrix::Matrix};
675    ///
676    /// let m = Matrix::identity() * 5.0;
677    /// let m_inv = m.inverse();
678    ///
679    /// assert_eq!(m_inv, Some(Matrix::with_diagonal([1.0 / 5.0; 3])));
680    /// ```
681    #[inline]
682    fn inverse(&self) -> Option<Self> {
683        #[inline]
684        fn cross(u: [f64; 3], v: [f64; 3]) -> [f64; 3] {
685            [
686                u[1] * v[2] - u[2] * v[1],
687                u[2] * v[0] - u[0] * v[2],
688                u[0] * v[1] - u[1] * v[0],
689            ]
690        }
691        let [x0, x1, x2] = self.rows;
692        let det = self.determinant();
693        if det == 0.0 {
694            return None;
695        }
696        let rows = [cross(x1, x2), cross(x2, x0), cross(x0, x1)];
697        Some(det.recip() * Self { rows }.transpose())
698    }
699}
700impl Invertible for Matrix<4, 4> {
701    /// Compute the inverse of a matrix. Will be `None` if the matrix is not invertible.
702    ///
703    /// This implementation uses a closed form solution for the matrix inverse based on
704    /// the Cayley–Hamilton method.
705    ///
706    /// # Examples
707    /// ```
708    /// use hoomd_linear_algebra::{
709    ///     Invertible, MatMul, SquareMatrix, matrix::Matrix44,
710    /// };
711    ///
712    /// let m = Matrix44::identity();
713    /// let m_inv = m.inverse().unwrap();
714    /// assert_eq!(m_inv.rows, m.rows);
715    /// ```
716    #[inline]
717    fn inverse(&self) -> Option<Self> {
718        let det = self.determinant();
719        if det == 0.0 {
720            return None;
721        }
722        // Compute components of Cayley–Hamilton factorization
723        let tr_a = self.trace();
724        let a_sq = self.powi(2);
725        let tr_a_sq = a_sq.trace();
726        let a_cb = a_sq.matmul(self);
727        let tr_a_cb = a_cb.trace();
728        let left =
729            (1.0 / 6.0) * (tr_a.powi(3) - 3.0 * tr_a * tr_a_sq + 2.0 * tr_a_cb) * Self::identity();
730        let center = (1.0 / 2.0) * *self * (tr_a.powi(2) - tr_a_sq);
731        Some(det.recip() * (left - center + a_sq * tr_a - a_cb))
732    }
733}
734
735impl<const N: usize, const M: usize> fmt::Display for Matrix<N, M> {
736    #[inline]
737    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
738        f.write_str(&format!(
739            "[{}]",
740            self.iter_rows()
741                .map(|row| {
742                    format!(
743                        "[{}]",
744                        row.iter()
745                            .map(ToString::to_string)
746                            .collect::<Vec<_>>()
747                            .join(", ")
748                    )
749                })
750                .collect::<Vec<_>>()
751                .join(",\n ")
752        ))
753    }
754}
755impl Matrix<2, 2> {
756    /// Decompose a [`Matrix22`] into a rotation, a scaling, and a second rotation.
757    ///
758    /// ```math
759    /// \mathbf{A} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^\intercal
760    /// ```
761    /// This implementation is based on the math in [Blinn 1996], and
762    /// ensures good (but not optimal) numerical stability. For certain
763    /// pathological inputs, preconditioning the matrix could provide a benefit
764    /// in numerical stability.
765    ///
766    /// `svd` sets all singular values to be positive.
767    ///
768    /// [Blinn 1996]: https://doi.org/10.1109/38.486688
769    ///
770    /// # Examples
771    /// ```
772    /// use hoomd_linear_algebra::{
773    ///     MatMul,
774    ///     matrix::{DiagonalMatrix, Matrix22},
775    /// };
776    /// let m = Matrix22 {
777    ///     rows: [[1.0, 2.0], [3.0, 4.0]],
778    /// };
779    /// let (u, s, vt) = m.svd();
780    /// let m_recon = u.matmul(&s.to_dense()).matmul(&vt);
781    /// for i in 0..2 {
782    ///     for j in 0..2 {
783    ///         assert!((m.rows[i][j] - m_recon.rows[i][j]).abs() < 1e-9);
784    ///     }
785    /// }
786    /// ```
787    #[must_use]
788    #[inline]
789    pub fn svd(&self) -> (Self, DiagonalMatrix<2>, Self) {
790        let a_plus_d = f64::midpoint(self[(0, 0)], self[(1, 1)]);
791
792        let a_minus_d = (self[(0, 0)] - self[(1, 1)]) / 2.0;
793        let b_plus_c = f64::midpoint(self[(0, 1)], self[(1, 0)]);
794        let b_minus_c = (self[(1, 0)] - self[(0, 1)]) / 2.0;
795        let (q, r) = (
796            (a_plus_d.powi(2) + b_minus_c.powi(2)).sqrt(),
797            (a_minus_d.powi(2) + b_plus_c.powi(2)).sqrt(),
798        );
799
800        let sy = q - r;
801        let sign_sy = sy.signum();
802
803        let (a1, a2) = (
804            f64::atan2(b_plus_c, a_minus_d),
805            f64::atan2(b_minus_c, a_plus_d),
806        );
807
808        let gamma = f64::midpoint(a1, a2);
809        let beta = (a2 - a1) / 2.0;
810
811        let (sr, cr) = beta.sin_cos();
812        let (sl, cl) = gamma.sin_cos();
813
814        let u = Matrix22 {
815            rows: [[cl, -sl], [sl, cl]],
816        };
817        let vt = Matrix22 {
818            rows: [[cr, -sr], [sr * sign_sy, cr * sign_sy]],
819        };
820
821        let singular_values = DiagonalMatrix::<2> {
822            elements: [q + r, sy.abs()],
823        };
824
825        (u, singular_values, vt)
826    }
827}
828
829impl Matrix<3, 3> {
830    /// Decompose a [`Matrix33`] into a rotation, a scaling, and a second rotation.
831    ///
832    /// ```math
833    /// \mathbf{A} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^\top
834    /// ```
835    /// This implementation is based on the method described by [McAdams 2011], which
836    ///  is a fast variant of the Jacobi iteration method.
837    ///
838    /// The method ensures that U and V are pure rotation matrices (determinant = 1).
839    /// As a result, the third singular value may be negative. For a conventional
840    /// SVD with non-negative singular values, the sign can be absorbed into U.
841    ///
842    /// [McAdams 2011]: https://digital.library.wisc.edu/1793/60736
843    ///
844    /// # Examples
845    /// ```
846    /// use hoomd_linear_algebra::{
847    ///     MatMul, SquareMatrix,
848    ///     matrix::{DiagonalMatrix, Matrix33},
849    /// };
850    /// let m = Matrix33 {
851    ///     rows: [[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]],
852    /// };
853    /// let (u, s, vt) = m.svd();
854    /// let m_recon = u.matmul(&s.to_dense()).matmul(&vt);
855    /// for i in 0..3 {
856    ///     for j in 0..3 {
857    ///         assert!((m.rows[i][j] - m_recon.rows[i][j]).abs() < 1e-9);
858    ///     }
859    /// }
860    /// ```
861    #[must_use]
862    #[inline]
863    pub fn svd(&self) -> (Self, DiagonalMatrix<3>, Self) {
864        #[inline]
865        fn jacobi_rotation(p_idx: usize, q_idx: usize, s: &mut Matrix33, v: &mut Matrix33) {
866            // ApproxGivensQuaternion adapted to build rotation matrix
867
868            let c_h = 2.0 * (s[(p_idx, p_idx)] - s[(q_idx, q_idx)]);
869            let s_h = s[(p_idx, q_idx)];
870            let gamma = 5.828_427_124_746_2; // 3 + 2 * sqrt(2)
871
872            let b = gamma * s_h.powi(2) < c_h.powi(2);
873            let omega = (c_h.powi(2) + s_h.powi(2)).sqrt().recip();
874
875            let (c_h_res, s_h_res) = if b {
876                (omega * c_h, omega * s_h)
877            } else {
878                (C_STAR, S_STAR)
879            };
880
881            // Build normalized rotation matrix from quaternion components
882            let c_h_sq = c_h_res.powi(2);
883            let s_h_sq = s_h_res.powi(2);
884            let scale = c_h_sq + s_h_sq;
885            let cos = (c_h_sq - s_h_sq) / scale;
886            let sin = (2.0 * c_h_res * s_h_res) / scale;
887
888            let mut q_mat = Matrix33::identity();
889            q_mat[(p_idx, p_idx)] = cos;
890            q_mat[(p_idx, q_idx)] = -sin;
891            q_mat[(q_idx, p_idx)] = sin;
892            q_mat[(q_idx, q_idx)] = cos;
893
894            *s = q_mat.transpose().matmul(s).matmul(&q_mat);
895            *v = v.matmul(&q_mat);
896        }
897
898        #[inline]
899        fn cond_neg_swap_rows(c: bool, rows: &mut [[f64; 3]], i: usize, j: usize) {
900            if c {
901                rows.swap(i, j);
902                rows[j].iter_mut().for_each(|x| *x *= -1.0);
903            }
904        }
905
906        #[inline]
907        fn qr_givens_rotation(p: usize, q: usize, r_mat: &mut Matrix33, u_mat: &mut Matrix33) {
908            let rho = (r_mat[(p, p)].powi(2) + r_mat[(q, p)].powi(2)).sqrt();
909            let c = if rho == 0.0 { 1.0 } else { r_mat[(p, p)] / rho };
910            let s_val = if rho == 0.0 { 0.0 } else { r_mat[(q, p)] / rho };
911
912            let mut q_t = Matrix33::identity();
913            q_t[(p, p)] = c;
914            q_t[(p, q)] = s_val;
915            q_t[(q, p)] = -s_val;
916            q_t[(q, q)] = c;
917
918            *r_mat = q_t.matmul(r_mat);
919            *u_mat = u_mat.matmul(&q_t.transpose());
920        }
921
922        // Symmetric Eigenanalysis of A^\top * A using Jacobi iterations
923        const NUM_JACOBI_SWEEPS: usize = 6; // Paper suggests 4, we want more accuracy
924        let mut singular_values = self.transpose().matmul(self);
925        let mut v = Self::identity();
926
927        for _ in 0..NUM_JACOBI_SWEEPS {
928            jacobi_rotation(0, 1, &mut singular_values, &mut v);
929            jacobi_rotation(0, 2, &mut singular_values, &mut v);
930            jacobi_rotation(1, 2, &mut singular_values, &mut v);
931        }
932
933        // Sort singular values and vectors
934        let mut b = self.matmul(&v);
935        let mut b_cols = b.transpose().rows;
936        let mut v_cols = v.transpose().rows;
937        let mut rhos: [f64; 3] = std::array::from_fn(|i| b_cols[i].iter().map(|&x| x * x).sum());
938
939        if rhos[0] < rhos[1] {
940            rhos.swap(0, 1);
941            cond_neg_swap_rows(true, &mut b_cols, 0, 1);
942            cond_neg_swap_rows(true, &mut v_cols, 0, 1);
943        }
944        if rhos[1] < rhos[2] {
945            rhos.swap(1, 2);
946            cond_neg_swap_rows(true, &mut b_cols, 1, 2);
947            cond_neg_swap_rows(true, &mut v_cols, 1, 2);
948        }
949        if rhos[0] < rhos[1] {
950            rhos.swap(0, 1);
951            cond_neg_swap_rows(true, &mut b_cols, 0, 1);
952            cond_neg_swap_rows(true, &mut v_cols, 0, 1);
953        }
954
955        b = Matrix { rows: b_cols }.transpose();
956        v = Matrix { rows: v_cols }.transpose();
957
958        // QR Decomposition of B = U * Sigma
959        let mut r = b;
960        let mut u = Self::identity();
961
962        qr_givens_rotation(0, 1, &mut r, &mut u);
963        qr_givens_rotation(0, 2, &mut r, &mut u);
964        qr_givens_rotation(1, 2, &mut r, &mut u);
965
966        // Enforce conventions for outputs. Rotations are proper and S can be negative
967        let mut sigma = r.diagonal();
968
969        if u.determinant() < 0.0 {
970            u.rows.iter_mut().for_each(|row| row[2] *= -1.0);
971            sigma[2] *= -1.0;
972        }
973
974        if v.determinant() < 0.0 {
975            v.rows.iter_mut().for_each(|row| row[2] *= -1.0);
976            sigma[2] *= -1.0;
977        }
978
979        (u, sigma, v.transpose())
980    }
981}
982
983#[cfg(test)]
984mod tests {
985    use std::{fmt::Debug, ops::Index};
986
987    use super::*;
988    use crate::matrix::{Matrix, Matrix22, Matrix33, Matrix44};
989    use approxim::{assert_relative_eq, assert_ulps_eq, ulps_eq};
990
991    use faer::Mat;
992    use rstest::rstest;
993
994    const EPS: f64 = 1e-13;
995
996    fn fill_faer<const N: usize, const M: usize>(m: [[f64; M]; N]) -> Mat<f64> {
997        let mut faer_matrix = Mat::<f64>::zeros(N, M);
998        for (i, row) in m.iter().enumerate() {
999            for (j, el) in row.iter().enumerate() {
1000                *faer_matrix.get_mut(i, j) = *el;
1001            }
1002        }
1003        faer_matrix
1004    }
1005    fn fill_faer_column<const N: usize>(c: [f64; N]) -> Mat<f64> {
1006        let mut faer_matrix = Mat::<f64>::zeros(N, 1);
1007        for (i, el) in c.iter().enumerate() {
1008            *faer_matrix.get_mut(i, 0) = *el;
1009        }
1010        faer_matrix
1011    }
1012    fn assert_matrices_ulps_eq<
1013        const N: usize,
1014        const M: usize,
1015        T0: Index<(usize, usize), Output = f64> + Debug,
1016        T1: Index<(usize, usize), Output = f64> + Debug,
1017    >(
1018        m0: &T0,
1019        m1: &T1,
1020    ) {
1021        for i in 0..N {
1022            for j in 0..M {
1023                if !ulps_eq!(m0[(i, j)], m1[(i, j)], epsilon = EPS) {
1024                    assert_ulps_eq!(m0[(i, j)], m1[(i, j)], epsilon = EPS);
1025                }
1026            }
1027        }
1028    }
1029    fn assert_diags_ulps_eq<const N: usize>(
1030        m0: &DiagonalMatrix<N>,
1031        m1: &impl std::ops::Index<usize, Output = f64>,
1032    ) {
1033        for i in 0..N {
1034            assert_ulps_eq!(m0[i], m1[i], epsilon = EPS);
1035        }
1036    }
1037    #[rstest(
1038        rows,
1039        case([[-9.0]]),
1040        case([[1.0, -2.0], [3.0, 4.0]]),
1041        case([[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]]),
1042        case([[2.0, 0.0, 1.0], [3.0, 9.0, 9.0], [5.0, 1.0, 1.0]]),
1043        case(Matrix::<4, 4>::identity().rows),
1044        case([
1045            [-10.0, 4.0, 3.0, 4.0],
1046            [300.0, 5.0, 6.0, 7.0],
1047            [3.0, 6.0, 8.0, 9.0],
1048            [4.0, 7.0, 9.0, 10.0]
1049        ]),
1050        case(Matrix::<5, 5>::full(3.6).diagonal().to_dense().rows),
1051        case(Matrix::<8, 8>::identity().rows),
1052    )]
1053    fn test_determinant<const N: usize>(rows: [[f64; N]; N]) {
1054        let matrix = Matrix { rows };
1055        let faer_matrix = fill_faer(rows);
1056
1057        let custom_det = matrix.determinant();
1058        let faer_det = faer_matrix.determinant();
1059
1060        assert_relative_eq!(custom_det, faer_det, max_relative = 1e-14);
1061    }
1062    #[rstest(
1063        a_rows, b_rows,
1064        case([[-9.0]], [[-9.0]]),
1065        case(
1066            [[1.0, -2.0], [3.0, 4.0]], [[0.0, 1.0], [1.0, 0.0]]
1067        ),
1068        case(
1069            [[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]],
1070            [[-2.0, 1.0, 0.0], [3.0, 0.0, 1.0], [1.0, 4.0, -1.0]]
1071        ),
1072        case(
1073            [[2.0, 0.0, 1.0], [3.0, 0.0, 0.0], [5.0, 1.0, 1.0]],
1074            [[1.0, 0.0, 2.0], [0.0, 1.0, 1.0], [4.0, 0.0, 0.0]]
1075        ),
1076        case(Matrix::<4, 4>::identity().rows, Matrix::<4, 4>::full(2.0).rows),
1077        case(Matrix::<5, 5>::full(3.6).diagonal().to_dense().rows, Matrix::<5, 5>::identity().rows),
1078        case(Matrix::<8, 8>::identity().rows, Matrix::<8, 8>::full(1.5).rows),
1079    )]
1080    fn test_matrix_multiply_square<const N: usize>(a_rows: [[f64; N]; N], b_rows: [[f64; N]; N]) {
1081        let a = Matrix { rows: a_rows };
1082        let b = Matrix { rows: b_rows };
1083
1084        let faer_a = fill_faer(a_rows);
1085        let faer_b = fill_faer(b_rows);
1086
1087        let custom_prod = a.matmul(&b);
1088        let faer_prod = faer_a * faer_b;
1089        assert_matrices_ulps_eq::<N, N, _, _>(&custom_prod, &faer_prod);
1090    }
1091
1092    #[rstest]
1093    #[case(
1094        [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
1095        [[7.0, 8.0], [9.0, 10.0], [11.0, 12.0]],
1096    )]
1097    #[case(
1098        [[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]],
1099        [[2.0, 3.0, 4.0], [5.0, 6.0, 7.0]],
1100    )]
1101    #[case(
1102        [[1.0, 2.0]],
1103        [[3.0], [4.0]],
1104    )]
1105    #[case(
1106        [[2.0, 0.0, 1.0]],
1107        [[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]],
1108    )]
1109    fn test_rectangular_matrix_multiply<const M: usize, const K: usize, const N: usize>(
1110        #[case] a_rows: [[f64; M]; N],
1111        #[case] b_rows: [[f64; K]; M],
1112    ) {
1113        let a = Matrix { rows: a_rows };
1114        let b = Matrix { rows: b_rows };
1115
1116        let faer_a = fill_faer(a_rows);
1117        let faer_b = fill_faer(b_rows);
1118
1119        let custom_prod = a.matmul(&b);
1120        let faer_prod = faer_a * faer_b;
1121        assert_matrices_ulps_eq::<N, K, _, _>(&custom_prod, &faer_prod);
1122    }
1123
1124    #[rstest(
1125        rows,
1126        case::identity(Matrix22::identity().rows),
1127        case::mixed_sign([[1.0, -2.0], [3.0, 4.0]]),
1128        case::det_zero([[12.0, 2.0], [4.0, 0.0]]),
1129        case::large_range([[1000.0, 0.0], [0.0, 1e-4]]),
1130        case::jordan_block([[1.0, 1.0], [0.0, 1.0]]),
1131        case::full_ones(Matrix22::full(1.0).rows),
1132        case::shear([[1.0, 2.0], [0.0, 1.0]]),
1133        case::nilpotent([[0.0, 1.0], [0.0, 0.0]]),
1134        case::scaling([[2.0, 0.0], [0.0, 3.0]]),
1135        // The fast algorithm does not compute the correct result for these degenerate
1136        // cases. test_svd_2x2_nalgebra verifies we reproduce the result expected for
1137        // this algorithm.
1138        // case::reflect([[0.0, -1.0], [1.0, 0.0]]),
1139        // case::negative_identity((Matrix22::identity()*-1.0).rows),
1140        // case::anti_diagonal([[0.0, 1.0], [1.0, 0.0]]),
1141        // case::singular([[1.0, 2.0], [2.0, 4.0]]),
1142    )]
1143    fn test_svd_2x2_faer(rows: [[f64; 2]; 2]) {
1144        let matrix = Matrix22 { rows };
1145        let (u, s, vt) = matrix.svd();
1146
1147        // Verify we can rebuild A from UΣVt
1148        assert_matrices_ulps_eq::<2, 2, _, _>(&u.matmul(&s.to_dense()).matmul(&vt), &matrix);
1149
1150        // Test against faer
1151        let faer = fill_faer(rows);
1152        let faersvd = faer.svd().unwrap();
1153        let (mut faeru, faers, mut faerv) =
1154            (faersvd.U().to_owned(), faersvd.S(), faersvd.V().to_owned());
1155
1156        if faeru.determinant().signum() != u.determinant().signum() {
1157            faeru[(0, 1)] *= -1.0;
1158            faeru[(1, 1)] *= -1.0;
1159        }
1160        if faerv.determinant().signum() != vt.determinant().signum() {
1161            faerv[(0, 1)] *= -1.0;
1162            faerv[(1, 1)] *= -1.0;
1163        }
1164
1165        assert_matrices_ulps_eq::<2, 2, _, _>(&u, &faeru);
1166        assert_diags_ulps_eq(&s, &faers);
1167        // Note that faer returns V, not Vt
1168        assert_matrices_ulps_eq::<2, 2, _, _>(&vt, &faerv.transpose());
1169    }
1170
1171    #[rstest(
1172        rows,
1173        case::identity(Matrix22::identity().rows),
1174        case::mixed_sign([[1.0, -2.0], [3.0, 4.0]]),
1175        case::det_zero([[12.0, 2.0], [4.0, 0.0]]),
1176        case::large_range([[1000.0, 0.0], [0.0, 1e-4]]),
1177        case::jordan_block([[1.0, 1.0], [0.0, 1.0]]),
1178        case::full_ones(Matrix22::full(1.0).rows),
1179        case::shear([[1.0, 2.0], [0.0, 1.0]]),
1180        case::nilpotent([[0.0, 1.0], [0.0, 0.0]]),
1181        case::scaling([[2.0, 0.0], [0.0, 3.0]]),
1182        case::reflect([[0.0, -1.0], [1.0, 0.0]]), // Numerical stability
1183        case::negative_identity((Matrix22::identity()*-1.0).rows),
1184        case::anti_diagonal([[0.0, 1.0], [1.0, 0.0]]),
1185        case::singular([[1.0, 2.0], [2.0, 4.0]]),
1186    )]
1187    fn test_svd_2x2_nalgebra(rows: [[f64; 2]; 2]) {
1188        let matrix = Matrix22 { rows };
1189        let (u, s, vt) = matrix.svd();
1190
1191        // Verify we can rebuild A from UΣVt
1192        assert_matrices_ulps_eq::<2, 2, _, _>(&u.matmul(&s.to_dense()).matmul(&vt), &matrix);
1193
1194        // Test against nalgebra
1195        let na = nalgebra::Matrix2::from(rows).transpose();
1196        let nasvd = na.svd(true, true);
1197        let (nau, nas, navt) = (nasvd.u.unwrap(), nasvd.singular_values, nasvd.v_t.unwrap());
1198
1199        assert_matrices_ulps_eq::<2, 2, _, _>(&u, &nau);
1200        assert_diags_ulps_eq::<2>(&s, &nas);
1201        assert_matrices_ulps_eq::<2, 2, _, _>(&vt, &navt);
1202    }
1203
1204    #[rstest(
1205        rows,
1206        case::identity(Matrix33::identity().rows),
1207        case::general([[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]]),
1208        case::symmetric([[1.0, 2.0, 3.0], [2.0, 5.0, 6.0], [3.0, 6.0, 9.0]]),
1209        case::near_singular([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.000_000_1]]),
1210        case::scaled((Matrix33::identity() * 10.0).rows),
1211        case::from_paper_repo([[0.433_807, 0.269_185, 0.543_034], [0.440_339, 0.443_024, 0.166_492], [0.793_913, 0.125_443, 0.730_333]]), // See https://github.com/ericjang/svd3
1212    )]
1213    fn test_svd_3x3_faer(rows: [[f64; 3]; 3]) {
1214        let matrix = Matrix33 { rows };
1215        let (u, s, vt) = matrix.svd();
1216
1217        // Verify reconstruction
1218        let m_recon = u.matmul(&s).matmul(&vt);
1219        assert_matrices_ulps_eq::<3, 3, _, _>(&m_recon, &matrix);
1220
1221        // Verify properties of U and V
1222        assert_relative_eq!(u.determinant(), 1.0, epsilon = EPS);
1223        assert_relative_eq!(vt.transpose().determinant(), 1.0, epsilon = EPS);
1224        assert_matrices_ulps_eq::<3, 3, _, _>(&u.matmul(&u.transpose()), &Matrix33::identity());
1225        assert_matrices_ulps_eq::<3, 3, _, _>(&vt.matmul(&vt.transpose()), &Matrix33::identity());
1226
1227        // Compare with faer SVD
1228        let faer_mat = fill_faer(rows);
1229        let faersvd = faer_mat.svd().unwrap();
1230
1231        let faers = faersvd.S();
1232        // Our implementation allows negative singular value
1233        assert_diags_ulps_eq(
1234            &DiagonalMatrix {
1235                elements: s.elements.map(f64::abs),
1236            },
1237            &faers,
1238        );
1239    }
1240
1241    #[test]
1242    fn test_transpose_2x2() {
1243        let rows = [[1.0, -2.0], [3.0, 4.0]];
1244        let matrix = Matrix::<2, 2> { rows };
1245        let faer_matrix = fill_faer(rows);
1246        let custom_transpose = matrix.transpose();
1247        let faer_transpose = faer_matrix.transpose();
1248        assert_matrices_ulps_eq::<2, 2, _, _>(&custom_transpose, &faer_transpose);
1249    }
1250
1251    #[test]
1252    fn test_transpose_2x3() {
1253        let rows = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
1254        let matrix = Matrix::<2, 3> { rows };
1255        let faer_matrix = fill_faer(rows);
1256        let custom_transpose = matrix.transpose();
1257        let faer_transpose = faer_matrix.transpose();
1258        assert_matrices_ulps_eq::<3, 2, _, _>(&custom_transpose, &faer_transpose);
1259    }
1260
1261    #[test]
1262    fn test_transpose_3x2() {
1263        let rows = [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]];
1264        let matrix = Matrix::<3, 2> { rows };
1265        let faer_matrix = fill_faer(rows);
1266        let custom_transpose = matrix.transpose();
1267        let faer_transpose = faer_matrix.transpose();
1268        assert_matrices_ulps_eq::<2, 3, _, _>(&custom_transpose, &faer_transpose);
1269    }
1270
1271    #[test]
1272    fn test_transpose_1x1() {
1273        let rows = [[-9.0]];
1274        let matrix = Matrix::<1, 1> { rows };
1275        assert_matrices_ulps_eq::<1, 1, _, _>(&matrix.transpose(), &matrix);
1276    }
1277
1278    #[test]
1279    fn test_general_matrix_methods() {
1280        // Matrix
1281        let zeros = Matrix::<2, 3>::zeros();
1282        let full = Matrix::<2, 3>::full(7.5);
1283        for i in 0..2 {
1284            for j in 0..3 {
1285                assert_eq!(zeros[(i, j)], 0.0);
1286                assert_eq!(full[(i, j)], 7.5);
1287            }
1288        }
1289    }
1290
1291    #[test]
1292    fn test_square_matrix_methods() {
1293        let identity = Matrix::<3, 3>::identity();
1294        let expected = Matrix::<3, 3> {
1295            rows: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
1296        };
1297        assert_matrices_ulps_eq::<3, 3, _, _>(&identity, &expected);
1298    }
1299
1300    #[test]
1301    fn test_diag_conversions() {
1302        let mat = Matrix::<3, 3> {
1303            rows: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]],
1304        };
1305        let diag = mat.diagonal();
1306        let expected_diag = DiagonalMatrix {
1307            elements: [1.0, 5.0, 9.0],
1308        };
1309        assert_diags_ulps_eq(&diag, &expected_diag);
1310
1311        let from_diag = diag.to_dense();
1312        let expected_from_diag = Matrix {
1313            rows: [[1.0, 0.0, 0.0], [0.0, 5.0, 0.0], [0.0, 0.0, 9.0]],
1314        };
1315        assert_matrices_ulps_eq::<3, 3, _, _>(&from_diag, &expected_from_diag);
1316    }
1317
1318    #[rstest(
1319        rows, vars,
1320        case(
1321            [[1.0, 2.0], [3.0, 4.0]],
1322            [0.5, 1.5]
1323        ),
1324        case(
1325            [[2.0, 0.0, 1.0], [3.0, 0.0, 0.0], [5.0, 1.0, 1.0]],
1326            [1.0, 2.0, 3.0]
1327        ),
1328        case(
1329            [[-33.0, 2.0, 0.0, 1.0], [3.0, -45.0, 0.0, 0.0], [5.0, 0.0, 1.0, 1.0], [0.0, 0.0, 0.0, 1.0]],
1330            [1.0, 2.0, 3.0, 4.0]
1331        ),
1332    )]
1333    fn test_quadratic_form<const N: usize>(rows: [[f64; N]; N], vars: [f64; N]) {
1334        let matrix = Matrix { rows };
1335        let result = matrix.compute_quadratic_form(&vars);
1336        assert_relative_eq!(
1337            result,
1338            (fill_faer_column(vars).transpose() * fill_faer(rows) * fill_faer_column(vars))[(0, 0)],
1339            max_relative = 1e-14
1340        );
1341    }
1342
1343    #[rstest(
1344        rows,
1345        case([[1.0, -2.0], [3.0, 4.0]]),
1346        case([[10.0, 0.0], [0.0, 0.1]]),
1347        case([[1.0, 1.0], [0.0, 1.0]]),
1348    )]
1349    fn test_inverse_2x2(rows: [[f64; 2]; 2]) {
1350        let matrix = Matrix22 { rows };
1351        let inv_matrix = matrix.inverse().expect("invertible");
1352        let product = matrix.matmul(&inv_matrix);
1353        let identity = Matrix22::identity();
1354
1355        assert_matrices_ulps_eq::<2, 2, _, _>(&product, &identity);
1356    }
1357    #[rstest(
1358        rows,
1359        case(Matrix33::identity().rows),
1360        case([[1.0, -3.0, 4.5], [3.0, 4.0,5.0], [8.0, -9.3, 10.0]]),
1361        case([[2.0, 1.0, 0.0], [0.0, 1.0, 2.0], [1.0, 0.0, 1.0]]),
1362        case([[5.0, -2.0, 3.0], [1.0, 0.0, 4.0], [-1.0, 2.0, 1.0]])
1363    )]
1364    fn test_inverse_3x3(rows: [[f64; 3]; 3]) {
1365        let matrix = Matrix33 { rows };
1366        let inv_matrix = matrix.inverse().expect("invertible");
1367        let product = matrix.matmul(&inv_matrix);
1368        let identity = Matrix33::identity();
1369
1370        assert_matrices_ulps_eq::<3, 3, _, _>(&product, &identity);
1371    }
1372    #[rstest(
1373        rows,
1374        case(Matrix44::identity().rows),
1375        case([[1.0, -4.0, 4.5,1.0], [4.0, 4.0,5.0,0.0], [8.0, -9.4, 10.0,9.0], [-1.0,-1.0,1.0,1.0]]),
1376    )]
1377    fn test_inverse_4x4(rows: [[f64; 4]; 4]) {
1378        let matrix = Matrix44 { rows };
1379        let inv_matrix = matrix.inverse().expect("invertible");
1380        let product = matrix.matmul(&inv_matrix);
1381        let identity = Matrix44::identity();
1382
1383        assert_matrices_ulps_eq::<4, 4, _, _>(&product, &identity);
1384    }
1385}
hoomd_linear_algebra/matrix.rs

hoomd_linear_algebra/
matrix.rs
