hoomd_linear_algebra/
lib.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
4//! Linear algebra optimized for small matrices.
5//!
6//! This crate places an emphasis on generality and simplicity, with optimization efforts
7//! targeted at small matrices. Some complex routines (SVD, matrix inversion, etc.) will
8//! only be implemented for certain shapes, and generally consist of specialized algorithms
9//! optimal for those inputs.
10//!
11//! This crate should not be considered a replacement for a dedicated linear
12//! algebra library like [faer-rs] or [nalgebra]. Instead, it can be used as a
13//! simple and lightweight dependency for matrix methods common to molecular
14//! simulation and analysis.
15//!
16//! [faer-rs]: https://github.com/sarah-quinones/faer-rs.git
17//! [nalgebra]: https://github.com/dimforge/nalgebra
18//!
19//! ## Matrices
20//!
21//! Construct a small rectangular [`Matrix`] on the stack:
22//! ```
23//! use hoomd_linear_algebra::matrix::Matrix;
24//!
25//! let m = Matrix {
26//!     rows: [[1.0, 2.0, 3.0], [-6.0, 3.0, 2.0]],
27//! };
28//! ```
29//!
30//! [`Matrix`]: matrix::Matrix
31//!
32//! Diagonal matrices ([`DiagonalMatrix`]) store only the diagonal elements:
33//! ```
34//! use hoomd_linear_algebra::matrix::DiagonalMatrix;
35//!
36//! let m = DiagonalMatrix {
37//!     elements: [-2.0, 4.0, -5.0],
38//! };
39//! ```
40//!
41//! [`DiagonalMatrix`]: matrix::DiagonalMatrix
42//!
43//! [`Matrix22`], [`Matrix33`], and [`Matrix44`] are type aliases for commonly used
44//! matrix sizes. Construct a 2x2 matrix with every element set to 4:
45//! ```
46//! use hoomd_linear_algebra::{Full, matrix::Matrix22};
47//!
48//! let m = Matrix22::full(4.0);
49//! ```
50//!
51//! [`Matrix22`]: matrix::Matrix22
52//! [`Matrix33`]: matrix::Matrix33
53//! [`Matrix44`]: matrix::Matrix44
54//!
55//! Construct a 3x3 identity matrix $` \mathbf{I} `$:
56//! ```
57//! use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix44};
58//!
59//! let m = Matrix44::identity();
60//! ```
61//!
62//! Index matrix entries by `(row, column)`:
63//! ```
64//! use hoomd_linear_algebra::matrix::Matrix;
65//!
66//! let m = Matrix {
67//!     rows: [[1.0, 2.0, 3.0], [-6.0, 3.0, 2.0]],
68//! };
69//!
70//! let element = m[(1, 2)];
71//! ```
72//!
73//! ## Matrix Operations
74//!
75//! Matrices can be added:
76//! ```
77//! use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
78//! let a = Matrix {
79//!     rows: [[1.0, -3.0], [-2.0, 4.0]],
80//! };
81//! let b = Matrix {
82//!     rows: [[4.0, -8.0], [6.0, 7.0]],
83//! };
84//!
85//! let mut c = a + b;
86//! c += a;
87//! ```
88//!
89//! subtracted:
90//! ```
91//! use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
92//! let a = Matrix {
93//!     rows: [[1.0, -3.0], [-2.0, 4.0]],
94//! };
95//! let b = Matrix {
96//!     rows: [[4.0, -8.0], [6.0, 7.0]],
97//! };
98//!
99//! let mut c = a - b;
100//! c += a;
101//! ```
102//!
103//! multiplied by a scalar:
104//! ```
105//! use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
106//! let a = Matrix {
107//!     rows: [[4.0, -8.0], [6.0, 7.0]],
108//! };
109//! let b = -2.0;
110//!
111//! let mut c = a * b;
112//! c *= b;
113//! ```
114//!
115//! negated:
116//! ```
117//! use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
118//! let a = Matrix {
119//!     rows: [[1.0, -3.0], [-2.0, 4.0]],
120//! };
121//!
122//! let b = -a;
123//! ```
124//!
125//! and indexed (in row,column ordering):
126//! ```
127//! use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
128//! let a = Matrix {
129//!     rows: [[1.0, -3.0], [-2.0, 4.0]],
130//! };
131//!
132//! let element = a[(1, 0)];
133//! ```
134//!
135//! You can also perform matrix-matrix multiplication:
136//! ```
137//! use hoomd_linear_algebra::{MatMul, matrix::Matrix};
138//!
139//! let a = Matrix {
140//!     rows: [[1.0, 2.0], [3.0, 4.0]],
141//! };
142//! let b = Matrix {
143//!     rows: [[4.0, 3.0], [2.0, 1.0]],
144//! };
145//!
146//! let c = a.matmul(&b);
147//! ```
148//!
149//! and invert matrices:
150//! ```
151//! use hoomd_linear_algebra::{Invertible, matrix::Matrix};
152//!
153//! let a = Matrix {
154//!     rows: [[1.0, 2.0], [3.0, 4.0]],
155//! };
156//!
157//! let b = a.inverse();
158//! ```
159//!
160//! ## Numerical Algorithms
161//!
162//! `hoomd-linear-algebra` implements a number of numerical algorithms on
163//! matrices:
164//!
165//! * [`Determinant`](matrix::Matrix::determinant)
166//! * [`Singular value decomposition (2x2)`](matrix::Matrix22::svd)
167//! * [`Singular value decomposition (3x3)`](matrix::Matrix33::svd)
168//! * [`Quadratic form`](QuadraticForm)
169//!
170//! # Complete documentation
171//!
172//! `hoomd-linear-algebra` is is a part of *hoomd-rs*. Read the [complete documentation]
173//! for more information.
174//!
175//! [complete documentation]: https://hoomd-rs.readthedocs.io
176
177use std::ops::{Add, AddAssign, Index, Mul, MulAssign, Neg, Sub, SubAssign};
178
179/// Structs implementing a large subset of Matrix traits.
180pub mod matrix;
181
182/// A lightweight representation of a diagonal matrix.
183mod diagonal;
184
185/// Compute the inverse of a matrix.
186///
187/// A matrix $`A`$ has an inverse $`A^{-1}`$ such that $`AA^{-1} = A^{-1}A = I`$.
188pub trait Invertible
189where
190    Self: Sized,
191{
192    /// Compute the inverse of a matrix.
193    ///
194    /// Returns `None` when the matrix is not invertible.
195    ///
196    /// # Example
197    ///
198    /// ```
199    /// use hoomd_linear_algebra::{Invertible, SquareMatrix, matrix::Matrix};
200    /// let m = Matrix::identity() * 5.0;
201    /// let m_inv = m.inverse();
202    ///
203    /// assert_eq!(m_inv, Some(Matrix::with_diagonal([1.0 / 5.0; 3])));
204    /// ```
205    #[must_use]
206    fn inverse(&self) -> Option<Self>;
207}
208
209/// Matrix multiplication.
210pub trait MatMul<Rhs> {
211    /// The type of the output matrix.
212    type Output;
213
214    /// Multiply two matrices.
215    ///
216    /// # Example
217    /// ```
218    /// use hoomd_linear_algebra::{
219    ///     Full, GeneralMatrix, MatMul, SquareMatrix,
220    ///     matrix::{DiagonalMatrix, Matrix, Matrix22},
221    /// };
222    ///
223    /// let a = Matrix22::full(5.0);
224    /// assert_eq!(a.matmul(&Matrix22::identity()), a);
225    ///
226    /// let b = Matrix::with_diagonal([3.0, 2.0]);
227    ///
228    /// assert_eq!(a.matmul(&b).rows, [[15.0, 10.0], [15.0, 10.0]]);
229    /// ```
230    #[must_use]
231    fn matmul(&self, rhs: &Rhs) -> Self::Output;
232}
233
234/// Common operations for all matrices.
235///
236/// Matrices can be added:
237/// ```
238/// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
239/// let a = Matrix {
240///     rows: [[1.0, -3.0], [-2.0, 4.0]],
241/// };
242/// let b = Matrix {
243///     rows: [[4.0, -8.0], [6.0, 7.0]],
244/// };
245///
246/// let mut c = a + b;
247/// c += a;
248/// ```
249///
250/// subtracted:
251/// ```
252/// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
253/// let a = Matrix {
254///     rows: [[1.0, -3.0], [-2.0, 4.0]],
255/// };
256/// let b = Matrix {
257///     rows: [[4.0, -8.0], [6.0, 7.0]],
258/// };
259///
260/// let mut c = a - b;
261/// c += a;
262/// ```
263///
264/// multiplied by a scalar:
265/// ```
266/// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
267/// let a = Matrix {
268///     rows: [[4.0, -8.0], [6.0, 7.0]],
269/// };
270/// let b = -2.0;
271///
272/// let mut c = a * b;
273/// c *= b;
274/// ```
275///
276/// negated:
277/// ```
278/// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
279/// let a = Matrix {
280///     rows: [[1.0, -3.0], [-2.0, 4.0]],
281/// };
282///
283/// let b = -a;
284/// ```
285///
286/// and indexed (in row,column ordering):
287/// ```
288/// use hoomd_linear_algebra::{Full, GeneralMatrix, matrix::Matrix};
289/// let a = Matrix {
290///     rows: [[1.0, -3.0], [-2.0, 4.0]],
291/// };
292///
293/// let element = a[(1, 0)];
294/// ```
295pub trait GeneralMatrix:
296    Add<Self, Output = Self>
297    + AddAssign<Self>
298    + Index<(usize, usize), Output = f64>
299    + Mul<f64, Output = Self>
300    + MulAssign<f64>
301    + Neg<Output = Self>
302    + Sized
303    + Sub<Self, Output = Self>
304    + SubAssign<Self>
305{
306    /// Fill a matrix with zeros.
307    ///
308    /// # Example
309    ///
310    /// ```
311    /// use hoomd_linear_algebra::{GeneralMatrix, matrix::Matrix};
312    ///
313    /// let a = Matrix::zeros();
314    ///
315    /// assert_eq!(a.rows, [[0.0, 0.0], [0.0, 0.0]]);
316    /// ```
317    #[must_use]
318    fn zeros() -> Self;
319
320    /// Get the shape of a matrix (``n_rows,n_columns``).
321    ///
322    /// # Example
323    ///
324    /// ```
325    /// use hoomd_linear_algebra::{GeneralMatrix, matrix::Matrix};
326    ///
327    /// let a = Matrix::<5, 7>::zeros();
328    /// assert_eq!(a.shape(), (5, 7));
329    /// ```
330    #[must_use]
331    fn shape(&self) -> (usize, usize);
332}
333
334/// Initialize matrices with identical elements.
335pub trait Full {
336    /// Construct a matrix the same value in every element.
337    ///
338    /// # Examples
339    /// ```
340    /// use hoomd_linear_algebra::{Full, matrix::Matrix22};
341    /// let m = Matrix22::full(5.0);
342    /// assert_eq!(m.rows, [[5.0, 5.0], [5.0, 5.0]]);
343    /// ```
344    #[must_use]
345    fn full(value: f64) -> Self;
346}
347
348/// Matrices that have the same number of rows and columns.
349pub trait SquareMatrix: GeneralMatrix {
350    /// Construct an N x N identity matrix.
351    ///
352    /// # Example
353    /// ```
354    /// use hoomd_linear_algebra::{SquareMatrix, matrix::Matrix22};
355    /// let m = Matrix22::identity();
356    /// assert_eq!(m.rows, [[1.0, 0.0], [0.0, 1.0]]);
357    /// ```
358    #[must_use]
359    fn identity() -> Self;
360}
361
362/// Solve the quadratic form.
363///
364/// ```math
365/// \mathbf{x}^{\intercal} \mathbf{A} \mathbf{x}
366/// ```
367pub trait QuadraticForm<const N: usize>: SquareMatrix {
368    /// Evaluate the quadratic form.
369    ///
370    /// The matrix `A` is given by `self` and the vector `x` in the argument.
371    #[must_use]
372    fn compute_quadratic_form(&self, x: &[f64; N]) -> f64;
373}
hoomd_linear_algebra/lib.rs

hoomd_linear_algebra/
lib.rs
