IMPL: Add trace, Hermitian conjugate, real/imag extractors (#87)

Add the small ergonomic operations requested in #87:

* Matrix::trace() and ComplexMatrix::trace(): sum of the diagonal,
  panics on non-square input (same contract as diag()).
* ComplexMatrix::h(): Hermitian conjugate (conjugate transpose).
* ComplexMatrix::real() / imag(): extract the real or imaginary part
  as a real Matrix, preserving the storage order.

The remaining items from #87 (IntoIterator support and a true matrix
exponential) are larger and stay tracked on the issue.

Author: Axect

src/complex/matrix.rs

@@ -12,7 +12,7 @@ use rand_distr::num_traits::{One, Zero};
 
 use crate::{
     complex::C64,
-    structure::matrix::Shape,
+    structure::matrix::{matrix, Matrix, Shape},
     traits::fp::{FPMatrix, FPVector},
     traits::general::Algorithm,
     traits::math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector},
@@ -74,6 +74,91 @@ impl ComplexMatrix {
     pub fn into_vec(self) -> Vec<C64> {
         self.data
     }
+
+    /// Trace of a square complex matrix (sum of the diagonal entries)
+    ///
+    /// # Panics
+    /// Panics if the matrix is not square.
+    ///
+    /// # Examples
+    /// ```rust
+    /// use peroxide::fuga::*;
+    /// use peroxide::complex::matrix::cmatrix;
+    ///
+    /// let a = cmatrix(vec![
+    ///     C64::new(1f64, 1f64),
+    ///     C64::new(2f64, 2f64),
+    ///     C64::new(3f64, 3f64),
+    ///     C64::new(4f64, -1f64),
+    /// ], 2, 2, Row);
+    /// assert_eq!(a.trace(), C64::new(5f64, 0f64));
+    /// ```
+    pub fn trace(&self) -> C64 {
+        self.diag().into_iter().sum()
+    }
+
+    /// Hermitian conjugate (conjugate transpose)
+    ///
+    /// # Examples
+    /// ```rust
+    /// use peroxide::fuga::*;
+    /// use peroxide::complex::matrix::cmatrix;
+    ///
+    /// let a = cmatrix(vec![
+    ///     C64::new(1f64, 1f64),
+    ///     C64::new(2f64, -1f64),
+    ///     C64::new(0f64, 2f64),
+    ///     C64::new(3f64, 0f64),
+    /// ], 2, 2, Row);
+    /// let ah = a.h();
+    /// assert_eq!(ah[(0, 1)], C64::new(0f64, -2f64));
+    /// assert_eq!(ah[(1, 0)], C64::new(2f64, 1f64));
+    /// ```
+    pub fn h(&self) -> Self {
+        let mut m = self.transpose();
+        for x in m.as_mut_slice() {
+            *x = x.conj();
+        }
+        m
+    }
+
+    /// Real part as a real matrix
+    ///
+    /// # Examples
+    /// ```rust
+    /// use peroxide::fuga::*;
+    /// use peroxide::complex::matrix::cmatrix;
+    ///
+    /// let a = cmatrix(vec![C64::new(1f64, 2f64), C64::new(3f64, 4f64)], 1, 2, Row);
+    /// assert_eq!(a.real()[(0, 1)], 3f64);
+    /// ```
+    pub fn real(&self) -> Matrix {
+        matrix(
+            self.data.iter().map(|c| c.re).collect::<Vec<f64>>(),
+            self.row,
+            self.col,
+            self.shape,
+        )
+    }
+
+    /// Imaginary part as a real matrix
+    ///
+    /// # Examples
+    /// ```rust
+    /// use peroxide::fuga::*;
+    /// use peroxide::complex::matrix::cmatrix;
+    ///
+    /// let a = cmatrix(vec![C64::new(1f64, 2f64), C64::new(3f64, 4f64)], 1, 2, Row);
+    /// assert_eq!(a.imag()[(0, 1)], 4f64);
+    /// ```
+    pub fn imag(&self) -> Matrix {
+        matrix(
+            self.data.iter().map(|c| c.im).collect::<Vec<f64>>(),
+            self.row,
+            self.col,
+            self.shape,
+        )
+    }
 }
 
 // =============================================================================

src/structure/matrix.rs

@@ -728,6 +728,22 @@ impl Matrix {
     pub fn into_vec(self) -> Vec<f64> {
         self.data
     }
+
+    /// Trace of a square matrix (sum of the diagonal entries)
+    ///
+    /// # Panics
+    /// Panics if the matrix is not square.
+    ///
+    /// # Examples
+    /// ```
+    /// use peroxide::fuga::*;
+    ///
+    /// let a = ml_matrix("1 2; 3 4");
+    /// assert_eq!(a.trace(), 5f64);
+    /// ```
+    pub fn trace(&self) -> f64 {
+        self.diag().iter().sum()
+    }
 }
 
 // =============================================================================

tests/complex_matrix.rs

@@ -21,3 +21,48 @@ fn test_seq() {
     );
     assert_eq!(v1.as_slice()[0], Complex64::new(1f64, 1f64));
 }
+
+#[cfg(feature = "complex")]
+#[test]
+fn test_trace_hermitian_real_imag() {
+    use num_complex::Complex64;
+    use peroxide::complex::matrix::cmatrix;
+
+    let a = cmatrix(
+        vec![
+            Complex64::new(1f64, 1f64),
+            Complex64::new(2f64, -1f64),
+            Complex64::new(0f64, 2f64),
+            Complex64::new(4f64, -1f64),
+        ],
+        2,
+        2,
+        Row,
+    );
+
+    // trace = (1+i) + (4-i) = 5
+    assert_eq!(a.trace(), Complex64::new(5f64, 0f64));
+
+    // hermitian conjugate: transpose + conjugate
+    let ah = a.h();
+    assert_eq!(ah[(0, 0)], Complex64::new(1f64, -1f64));
+    assert_eq!(ah[(0, 1)], Complex64::new(0f64, -2f64));
+    assert_eq!(ah[(1, 0)], Complex64::new(2f64, 1f64));
+    assert_eq!(ah[(1, 1)], Complex64::new(4f64, 1f64));
+
+    // (A^H)^H == A
+    let ahh = ah.h();
+    for i in 0..2 {
+        for j in 0..2 {
+            assert_eq!(ahh[(i, j)], a[(i, j)]);
+        }
+    }
+
+    // real / imag parts
+    let re = a.real();
+    let im = a.imag();
+    assert_eq!(re[(0, 1)], 2f64);
+    assert_eq!(im[(0, 1)], -1f64);
+    assert_eq!(re[(1, 0)], 0f64);
+    assert_eq!(im[(1, 0)], 2f64);
+}

tests/matrix.rs

@@ -170,3 +170,11 @@ fn test_matrix_rejects_length_mismatch() {
     // row * col must be impossible to construct from safe code.
     let _ = matrix(vec![1, 2, 3], 2, 2, Row);
 }
+
+#[test]
+fn test_trace() {
+    let a = ml_matrix("1 2; 3 4");
+    assert_eq!(a.trace(), 5f64);
+    // layout independence
+    assert_eq!(a.change_shape().trace(), 5f64);
+}