feat: Add Bunch–Kaufman LDL Decomposition

src/base/matrix_view.rs

@@ -33,14 +33,6 @@ macro_rules! view_storage_impl (
         #[deprecated = "Use ViewStorage(Mut) instead."]
         pub type $legacy_name<'a, T, R, C, RStride, CStride> = $T<'a, T, R, C, RStride, CStride>;
 
-        unsafe impl<'a, T: Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send
-            for $T<'a, T, R, C, RStride, CStride>
-        {}
-
-        unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync
-            for $T<'a, T, R, C, RStride, CStride>
-        {}
-
         impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> $T<'a, T, R, C, RStride, CStride> {
             /// Create a new matrix view without bounds checking and from a raw pointer.
             ///
@@ -136,6 +128,20 @@ impl<T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Copy
 {
 }
 
+/// Safety: Equivalent to a shared reference to `T`. All `Dim` type arguments are `Send + Sync`. A
+/// shared reference can be sent iff `T: Sync`.
+unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send
+    for ViewStorage<'a, T, R, C, RStride, CStride>
+{
+}
+
+/// Safety: Equivalent to a shared reference to `T`. All `Dim` type arguments are `Send + Sync`. A
+/// shared reference is `Sync` iff `T: Sync`.
+unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync
+    for ViewStorage<'a, T, R, C, RStride, CStride>
+{
+}
+
 impl<T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Clone
     for ViewStorage<'_, T, R, C, RStride, CStride>
 {
@@ -145,6 +151,20 @@ impl<T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Clone
     }
 }
 
+/// Safety: Equivalent to a unique reference to `T`. All `Dim` type arguments are `Send + Sync`. A
+/// unique reference is `Send` iff `T: Send`.
+unsafe impl<'a, T: Send, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Send
+    for ViewStorageMut<'a, T, R, C, RStride, CStride>
+{
+}
+
+/// Safety: Equivalent to a unique reference to `T`. All `Dim` type arguments are `Send + Sync`. A
+/// unique reference is `Sync` iff `T: Sync`.
+unsafe impl<'a, T: Sync, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Sync
+    for ViewStorageMut<'a, T, R, C, RStride, CStride>
+{
+}
+
 impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim>
     ViewStorageMut<'a, T, R, C, RStride, CStride>
 where

src/linalg/decomposition.rs

@@ -1,8 +1,8 @@
 use crate::storage::Storage;
 use crate::{
     Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff,
-    DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, LU, Matrix, OMatrix, QR, RealField, SVD,
-    Schur, SymmetricEigen, SymmetricTridiagonal, U1, UDU,
+    DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, LDL, LU, Matrix, OMatrix, QR, RealField,
+    SVD, Schur, SymmetricEigen, SymmetricTridiagonal, U1, UDU,
 };
 
 /// # Rectangular matrix decomposition
@@ -244,6 +244,7 @@ impl<T: ComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 /// | -------------------------|---------------------------|--------------|
 /// | Hessenberg               | `Q * H * Qᵀ`             | `Q` is a unitary matrix and `H` an upper-Hessenberg matrix. |
 /// | Cholesky                 | `L * Lᵀ`                 | `L` is a lower-triangular matrix. |
+/// | LDL decomposition        | `Pᵀ * L * D * Lᴴ * P`    | `L` is unit lower-triangular, `D` is Hermitian block-diagonal, and `P` is a permutation matrix. |
 /// | UDU                      | `U * D * Uᵀ`             | `U` is a upper-triangular matrix, and `D` a diagonal matrix. |
 /// | Schur decomposition      | `Q * T * Qᵀ`             | `Q` is an unitary matrix and `T` a quasi-upper-triangular matrix. |
 /// | Symmetric eigendecomposition | `Q ~ Λ ~ Qᵀ`   | `Q` is an unitary matrix, and `Λ` is a real diagonal matrix. |
@@ -260,6 +261,18 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> Matrix<T, D, D, S> {
         Cholesky::new(self.into_owned())
     }
 
+    /// Computes the LDL decomposition of this matrix.
+    /// The input matrix `self` is assumed to be Hermitian (symmetric) and the decomposition will
+    /// only read the lower-triangular part of `self`.
+    pub fn ldl(self) -> LDL<T, D>
+    where
+        T: Copy,
+        T::RealField: Copy,
+        DefaultAllocator: Allocator<D> + Allocator<D, D>,
+    {
+        LDL::new(self.into_owned())
+    }
+
     /// Attempts to compute the UDU decomposition of this matrix.
     ///
     /// The input matrix `self` is assumed to be symmetric and this decomposition will only read