panic docs and rank-deficiency robustness for fullpivlu and colpivqr

Cargo.toml

@@ -18,7 +18,7 @@ assert2 = "0.3"
 num-traits = "0.2"
 num-complex = "0.4.3"
 rayon = "1.7"
-pulp = { version = "0.13", default-features = false }
+pulp = { version = "0.13.2", default-features = false }
 bytemuck = "1"
 
 [profile.dev]

faer-cholesky/src/ldlt_diagonal/compute.rs

@@ -291,14 +291,17 @@ fn cholesky_in_place_impl<E: ComplexField>(
 ///
 /// # Warning
 ///
-/// The Cholesky decomposition may have poor numerical stability properties when used with non
-/// positive definite matrices. In the general case, it is recommended to first permute the matrix
-/// using [`crate::compute_cholesky_permutation`] and
-/// [`permute_rows_and_cols_symmetric`](faer_core::permutation::permute_rows_and_cols_symmetric_lower).
+/// The Cholesky decomposition with diagonal may have poor numerical stability properties when used
+/// with non positive definite matrices. In the general case, it is recommended to first permute
+/// (and conjugate when necessary) the rows and columns of the matrix using the permutation obtained
+/// from [`crate::compute_cholesky_permutation`].
 ///
 /// # Panics
 ///
 /// Panics if the input matrix is not square.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`raw_cholesky_in_place_req`]).
 #[track_caller]
 #[inline]
 pub fn raw_cholesky_in_place<E: ComplexField>(

faer-cholesky/src/ldlt_diagonal/mod.rs

@@ -3,9 +3,9 @@
 //! where $D$ is a diagonal matrix, and $L$ is a unit lower triangular matrix.
 //!
 //! The Cholesky decomposition with diagonal may have poor numerical stability properties when used
-//! with non positive definite matrices. In the general case, it is recommended to first permute the
-//! matrix using [`crate::compute_cholesky_permutation`] and
-//! [`permute_rows_and_cols_symmetric`](faer_core::permutation::permute_rows_and_cols_symmetric_lower).
+//! with non positive definite matrices. In the general case, it is recommended to first permute
+//! (and conjugate when necessary) the rows and columns of the matrix using the permutation obtained
+//! from [`crate::compute_cholesky_permutation`].
 
 pub mod compute;
 pub mod solve;
@@ -239,8 +239,9 @@ mod tests {
             delete_rows_and_cols_clobber(
                 a.as_mut(),
                 &mut [1, 3],
+                Parallelism::None,
                 DynStack::new(&mut GlobalMemBuffer::new(
-                    delete_rows_and_cols_clobber_req::<E>(n, r).unwrap(),
+                    delete_rows_and_cols_clobber_req::<E>(n, r, Parallelism::None).unwrap(),
                 )),
             );
 
@@ -268,8 +269,9 @@ mod tests {
             delete_rows_and_cols_clobber(
                 a.as_mut(),
                 &mut [0, 2],
+                Parallelism::None,
                 DynStack::new(&mut GlobalMemBuffer::new(
-                    delete_rows_and_cols_clobber_req::<E>(n, r).unwrap(),
+                    delete_rows_and_cols_clobber_req::<E>(n, r, Parallelism::None).unwrap(),
                 )),
             );
 
@@ -297,8 +299,9 @@ mod tests {
             delete_rows_and_cols_clobber(
                 a.as_mut(),
                 &mut [0, 2, 3],
+                Parallelism::None,
                 DynStack::new(&mut GlobalMemBuffer::new(
-                    delete_rows_and_cols_clobber_req::<E>(n, r).unwrap(),
+                    delete_rows_and_cols_clobber_req::<E>(n, r, Parallelism::None).unwrap(),
                 )),
             );
 

faer-cholesky/src/ldlt_diagonal/solve.rs

@@ -1,16 +1,78 @@
-use dyn_stack::DynStack;
-use faer_core::{solve, zipped, ComplexField, Conj, MatMut, MatRef, Parallelism};
+use dyn_stack::{DynStack, SizeOverflow, StackReq};
+use faer_core::{solve, zipped, ComplexField, Conj, Entity, MatMut, MatRef, Parallelism};
 use reborrow::*;
 
 use assert2::assert;
 
+/// Computes the size and alignment of required workspace for solving a linear system defined by a
+/// matrix in place, given its diagonal LDLT decomposition.
+pub fn solve_in_place_req<E: Entity>(
+    cholesky_dimension: usize,
+    rhs_ncols: usize,
+    parallelism: Parallelism,
+) -> Result<StackReq, SizeOverflow> {
+    let _ = cholesky_dimension;
+    let _ = rhs_ncols;
+    let _ = parallelism;
+    Ok(StackReq::default())
+}
+
+/// Computes the size and alignment of required workspace for solving a linear system defined by a
+/// matrix out of place, given its diagonal LDLT decomposition.
+pub fn solve_req<E: Entity>(
+    cholesky_dimension: usize,
+    rhs_ncols: usize,
+    parallelism: Parallelism,
+) -> Result<StackReq, SizeOverflow> {
+    let _ = cholesky_dimension;
+    let _ = rhs_ncols;
+    let _ = parallelism;
+    Ok(StackReq::default())
+}
+
+/// Computes the size and alignment of required workspace for solving a linear system defined by
+/// the transpose of a matrix in place, given its diagonal LDLT decomposition.
+pub fn solve_transpose_in_place_req<E: Entity>(
+    cholesky_dimension: usize,
+    rhs_ncols: usize,
+    parallelism: Parallelism,
+) -> Result<StackReq, SizeOverflow> {
+    let _ = cholesky_dimension;
+    let _ = rhs_ncols;
+    let _ = parallelism;
+    Ok(StackReq::default())
+}
+
+/// Computes the size and alignment of required workspace for solving a linear system defined by
+/// the transpose of a matrix out of place, given its diagonal LDLT decomposition.
+pub fn solve_transpose_req<E: Entity>(
+    cholesky_dimension: usize,
+    rhs_ncols: usize,
+    parallelism: Parallelism,
+) -> Result<StackReq, SizeOverflow> {
+    let _ = cholesky_dimension;
+    let _ = rhs_ncols;
+    let _ = parallelism;
+    Ok(StackReq::default())
+}
+
 /// Given the Cholesky factors of a matrix $A$ and a matrix $B$ stored in `rhs`, this function
 /// computes the solution of the linear system:
 /// $$\text{Op}_A(A)X = B.$$
 ///
 /// $\text{Op}_A$ is either the identity or the conjugation depending on the value of `conj_lhs`.  
 ///
 /// The solution of the linear system is stored in `rhs`.
+///
+/// # Panics
+///
+/// Panics if any of these conditions is violated:
+///
+/// * `cholesky_factors` must be square of dimension `n`.
+/// * `rhs` must have `n` rows.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`solve_in_place_req`]).
 #[track_caller]
 pub fn solve_in_place_with_conj<E: ComplexField>(
     cholesky_factors: MatRef<'_, E>,
@@ -60,6 +122,16 @@ pub fn solve_in_place_with_conj<E: ComplexField>(
 /// $\text{Op}_A$ is either the identity or the conjugation depending on the value of `conj_lhs`.  
 ///
 /// The solution of the linear system is stored in `rhs`.
+///
+/// # Panics
+///
+/// Panics if any of these conditions is violated:
+///
+/// * `cholesky_factors` must be square of dimension `n`.
+/// * `rhs` must have `n` rows.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`solve_transpose_in_place_req`]).
 #[track_caller]
 pub fn solve_transpose_in_place_with_conj<E: ComplexField>(
     cholesky_factors: MatRef<'_, E>,
@@ -88,6 +160,17 @@ pub fn solve_transpose_in_place_with_conj<E: ComplexField>(
 /// $\text{Op}_A$ is either the identity or the conjugation depending on the value of `conj_lhs`.  
 ///
 /// The solution of the linear system is stored in `dst`.
+///
+/// # Panics
+///
+/// Panics if any of these conditions is violated:
+///
+/// * `cholesky_factors` must be square of dimension `n`.
+/// * `rhs` must have `n` rows.
+/// * `dst` must have `n` rows, and the same number of columns as `rhs`.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`solve_transpose_req`]).
 #[track_caller]
 pub fn solve_transpose_with_conj<E: ComplexField>(
     dst: MatMut<'_, E>,
@@ -109,6 +192,17 @@ pub fn solve_transpose_with_conj<E: ComplexField>(
 /// $\text{Op}_A$ is either the identity or the conjugation depending on the value of `conj_lhs`.  
 ///
 /// The solution of the linear system is stored in `dst`.
+///
+/// # Panics
+///
+/// Panics if any of these conditions is violated:
+///
+/// * `cholesky_factors` must be square of dimension `n`.
+/// * `rhs` must have `n` rows.
+/// * `dst` must have `n` rows, and the same number of columns as `rhs`.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`solve_req`]).
 #[track_caller]
 pub fn solve_with_conj<E: ComplexField>(
     dst: MatMut<'_, E>,

faer-cholesky/src/ldlt_diagonal/update.rs

@@ -637,6 +637,14 @@ impl<'a, E: ComplexField> RankRUpdate<'a, E> {
 ///
 /// The matrix $W$ and the vector $\alpha$ are clobbered, meaning that the values they contain after
 /// the function is called are unspecified.
+///
+/// # Panics
+///
+/// Panics if any of these conditions is violated:
+/// * `cholesky_factors` must be square of dimension `n`.
+/// * `w` must have `n` rows.
+/// * `alpha` must have one column.
+/// * `alpha` must have the same number of rows as the number of columns in `w`.
 #[track_caller]
 pub fn rank_r_update_clobber<E: ComplexField>(
     cholesky_factors: MatMut<'_, E>,
@@ -726,7 +734,9 @@ pub(crate) fn rank_update_indices(
 pub fn delete_rows_and_cols_clobber_req<E: Entity>(
     dim: usize,
     number_of_rows_to_remove: usize,
+    parallelism: Parallelism,
 ) -> Result<StackReq, SizeOverflow> {
+    let _ = parallelism;
     let r = number_of_rows_to_remove;
     StackReq::try_all_of([temp_mat_req::<E>(dim, r)?, temp_mat_req::<E>(r, 1)?])
 }
@@ -740,12 +750,22 @@ pub fn delete_rows_and_cols_clobber_req<E: Entity>(
 ///
 /// The indices are clobbered, meaning that the values that the slice contains after the function
 /// is called are unspecified.
+///
+/// # Panics
+///
+/// Panics if the matrix is not square, the indices are out of bounds, or the list of indices
+/// contains duplicates.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`delete_rows_and_cols_clobber_req`]).
 #[track_caller]
 pub fn delete_rows_and_cols_clobber<E: ComplexField>(
     cholesky_factors: MatMut<'_, E>,
     indices: &mut [usize],
+    parallelism: Parallelism,
     stack: DynStack<'_>,
 ) {
+    let _ = parallelism;
     let n = cholesky_factors.nrows();
     let r = indices.len();
     assert!(cholesky_factors.ncols() == n);
@@ -824,6 +844,13 @@ pub fn insert_rows_and_cols_clobber_req<E: Entity>(
 ///
 /// The inserted matrix is clobbered, meaning that the values it contains after the function
 /// is called are unspecified.
+///
+/// # Panics
+///
+/// Panics if the index is out of bounds or the matrix dimensions are invalid.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`insert_rows_and_cols_clobber_req`]).
 #[track_caller]
 pub fn insert_rows_and_cols_clobber<E: ComplexField>(
     cholesky_factors_extended: MatMut<'_, E>,

faer-cholesky/src/llt/compute.rs

@@ -185,8 +185,10 @@ fn cholesky_in_place_impl<E: ComplexField>(
 ///
 /// # Panics
 ///
-/// - Panics if the input matrix is not square.
-/// - Panics if the provided memory in `stack` is insufficient.
+/// Panics if the input matrix is not square.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`cholesky_in_place_req`]).
 #[track_caller]
 #[inline]
 pub fn cholesky_in_place<E: ComplexField>(

faer-cholesky/src/llt/inverse.rs

@@ -65,7 +65,9 @@ pub fn invert_lower_in_place_req<E: Entity>(
 /// # Panics
 ///
 /// - Panics if `cholesky_factor` is not a square matrix.
-/// - Panics if the provided memory in `stack` is insufficient.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`invert_lower_req`]).
 #[track_caller]
 pub fn invert_lower_in_place<E: ComplexField>(
     cholesky_factor: MatMut<'_, E>,
@@ -83,7 +85,9 @@ pub fn invert_lower_in_place<E: ComplexField>(
 ///
 /// - Panics if `cholesky_factor` is not a square matrix.
 /// - Panics if the destination shape doesn't match the shape of the matrix.
-/// - Panics if the provided memory in `stack` is insufficient.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`invert_lower_in_place_req`]).
 #[track_caller]
 pub fn invert_lower<E: ComplexField>(
     dst: MatMut<'_, E>,

faer-cholesky/src/llt/mod.rs

@@ -243,8 +243,9 @@ mod tests {
             delete_rows_and_cols_clobber(
                 a.as_mut(),
                 &mut [1, 3],
+                Parallelism::Rayon(8),
                 DynStack::new(&mut GlobalMemBuffer::new(
-                    delete_rows_and_cols_clobber_req::<E>(n, r).unwrap(),
+                    delete_rows_and_cols_clobber_req::<E>(n, r, Parallelism::Rayon(8)).unwrap(),
                 )),
             );
 
@@ -270,8 +271,9 @@ mod tests {
             delete_rows_and_cols_clobber(
                 a.as_mut(),
                 &mut [0, 2],
+                Parallelism::Rayon(8),
                 DynStack::new(&mut GlobalMemBuffer::new(
-                    delete_rows_and_cols_clobber_req::<E>(n, r).unwrap(),
+                    delete_rows_and_cols_clobber_req::<E>(n, r, Parallelism::Rayon(8)).unwrap(),
                 )),
             );
 
@@ -297,8 +299,9 @@ mod tests {
             delete_rows_and_cols_clobber(
                 a.as_mut(),
                 &mut [0, 2, 3],
+                Parallelism::Rayon(8),
                 DynStack::new(&mut GlobalMemBuffer::new(
-                    delete_rows_and_cols_clobber_req::<E>(n, r).unwrap(),
+                    delete_rows_and_cols_clobber_req::<E>(n, r, Parallelism::Rayon(8)).unwrap(),
                 )),
             );
 
@@ -335,8 +338,10 @@ mod tests {
             delete_rows_and_cols_clobber(
                 a.as_mut(),
                 &mut [position, position + 1],
+                Parallelism::Rayon(8),
                 DynStack::new(&mut GlobalMemBuffer::new(
-                    delete_rows_and_cols_clobber_req::<f64>(n + r, r).unwrap(),
+                    delete_rows_and_cols_clobber_req::<f64>(n + r, r, Parallelism::Rayon(8))
+                        .unwrap(),
                 )),
             );
 

faer-cholesky/src/llt/reconstruct.rs

@@ -30,6 +30,9 @@ pub fn reconstruct_lower_in_place_req<E: Entity>(
 ///
 /// - Panics if `cholesky_factor` is not a square matrix.
 /// - Panics if the destination shape doesn't match the shape of the matrix.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`reconstruct_lower_req`]).
 #[track_caller]
 pub fn reconstruct_lower<E: ComplexField>(
     dst: MatMut<'_, E>,
@@ -60,6 +63,9 @@ pub fn reconstruct_lower<E: ComplexField>(
 /// # Panics
 ///
 /// - Panics if `cholesky_factor` is not a square matrix.
+///
+/// This can also panic if the provided memory in `stack` is insufficient (see
+/// [`reconstruct_lower_in_place_req`]).
 #[track_caller]
 pub fn reconstruct_lower_in_place<E: ComplexField>(
     cholesky_factor: MatMut<'_, E>,