Skip to main content

a/math/matrix/
basic.rs

1//! Transposes, products elementwise, row and column operations, the trace and powers.
2
3use super::{Element, Matrix};
4
5impl<T: Clone> Matrix<T> {
6    /// The transpose `Aᵀ`, whose element in row `i` and column `j` is this one's in row `j` and
7    /// column `i`.
8    pub fn transpose(&self) -> Self {
9        let rows = self.rows;
10        Matrix::from_fn(self.cols, rows, |i, j| self.data[j + i * rows].clone())
11    }
12}
13
14impl<T: Element> Matrix<T> {
15    /// The conjugate transpose `A*`, the transpose with each element conjugated: the transpose
16    /// itself for real elements.
17    pub fn conjugate_transpose(&self) -> Self {
18        let rows = self.rows;
19        Matrix::from_fn(self.cols, rows, |i, j| self.data[j + i * rows].conj())
20    }
21
22    /// The product elementwise, Hadamard's `A ∘ B`.
23    ///
24    /// # Panics
25    ///
26    /// Panics if the matrices are not of the same shape.
27    #[track_caller]
28    pub fn hadamard(&self, other: &Matrix<T>) -> Matrix<T> {
29        assert!(
30            self.shape() == other.shape(),
31            "the matrices are not of the same shape"
32        );
33        Matrix {
34            rows: self.rows,
35            cols: self.cols,
36            data: self
37                .data
38                .iter()
39                .zip(&other.data)
40                .map(|(a, b)| T::mul(a, b))
41                .collect(),
42        }
43    }
44
45    /// Swaps rows `a` and `b`.
46    ///
47    /// # Panics
48    ///
49    /// Panics if either is outside the matrix.
50    #[track_caller]
51    pub fn swap_rows(&mut self, a: usize, b: usize) {
52        assert!(
53            a < self.rows && b < self.rows,
54            "the row is outside the matrix"
55        );
56        if a != b {
57            for j in 0..self.cols {
58                self.data.swap(a + j * self.rows, b + j * self.rows);
59            }
60        }
61    }
62
63    /// Multiplies row `i` by `k`.
64    ///
65    /// # Panics
66    ///
67    /// Panics if the row is outside the matrix.
68    #[track_caller]
69    pub fn scale_row(&mut self, i: usize, k: &T) {
70        assert!(i < self.rows, "the row is outside the matrix");
71        for j in 0..self.cols {
72            T::mul_assign(&mut self.data[i + j * self.rows], k);
73        }
74    }
75
76    /// Adds `k` times row `source` to row `target`.
77    ///
78    /// # Panics
79    ///
80    /// Panics if either row is outside the matrix.
81    #[track_caller]
82    pub fn add_scaled_row(&mut self, target: usize, source: usize, k: &T) {
83        assert!(
84            target < self.rows && source < self.rows,
85            "the row is outside the matrix"
86        );
87        for j in 0..self.cols {
88            let x = T::mul(&self.data[source + j * self.rows], k);
89            T::add_assign(&mut self.data[target + j * self.rows], &x);
90        }
91    }
92
93    /// Swaps columns `a` and `b`.
94    ///
95    /// # Panics
96    ///
97    /// Panics if either is outside the matrix.
98    #[track_caller]
99    pub fn swap_columns(&mut self, a: usize, b: usize) {
100        assert!(
101            a < self.cols && b < self.cols,
102            "the column is outside the matrix"
103        );
104        if a != b {
105            let (lo, hi) = (a.min(b), a.max(b));
106            let (left, right) = self.data.split_at_mut(hi * self.rows);
107            left[lo * self.rows..(lo + 1) * self.rows].swap_with_slice(&mut right[..self.rows]);
108        }
109    }
110
111    /// Multiplies column `j` by `k`.
112    ///
113    /// # Panics
114    ///
115    /// Panics if the column is outside the matrix.
116    #[track_caller]
117    pub fn scale_column(&mut self, j: usize, k: &T) {
118        for x in self.column_mut(j) {
119            T::mul_assign(x, k);
120        }
121    }
122
123    /// Adds `k` times column `source` to column `target`.
124    ///
125    /// # Panics
126    ///
127    /// Panics if either column is outside the matrix.
128    #[track_caller]
129    pub fn add_scaled_column(&mut self, target: usize, source: usize, k: &T) {
130        assert!(
131            target < self.cols && source < self.cols,
132            "the column is outside the matrix"
133        );
134        let rows = self.rows;
135        for i in 0..rows {
136            let x = T::mul(&self.data[i + source * rows], k);
137            T::add_assign(&mut self.data[i + target * rows], &x);
138        }
139    }
140
141    /// The trace `tr(A)`, the sum of the main diagonal.
142    ///
143    /// # Panics
144    ///
145    /// Panics if the matrix is not square.
146    #[track_caller]
147    pub fn trace(&self) -> T {
148        assert!(self.is_square(), "the matrix is not square");
149        let mut t = T::zero();
150        for i in 0..self.rows {
151            T::add_assign(&mut t, &self.data[i + i * self.rows]);
152        }
153        t
154    }
155
156    /// `Aⁿ`, by repeated squaring: about `2 log₂ n` products. `A⁰` is the identity.
157    ///
158    /// # Panics
159    ///
160    /// Panics if the matrix is not square.
161    #[track_caller]
162    pub fn pow(&self, n: u32) -> Self {
163        assert!(self.is_square(), "the matrix is not square");
164        if n == 0 {
165            return Matrix::identity(self.rows);
166        }
167        let mut base = self.clone();
168        let mut n = n;
169        while n & 1 == 0 {
170            base = &base * &base;
171            n >>= 1;
172        }
173        let mut r = base.clone();
174        n >>= 1;
175        while n > 0 {
176            base = &base * &base;
177            if n & 1 == 1 {
178                r = &r * &base;
179            }
180            n >>= 1;
181        }
182        r
183    }
184}