Cholesky factorization of PSD matrices

src/cholesky.rs

@@ -48,6 +48,7 @@ use num_traits::Float;
 
 use crate::convert::*;
 use crate::error::*;
+use crate::lapack::cholesky_semi::cholesky_semi;
 use crate::layout::*;
 use crate::triangular::IntoTriangular;
 use crate::types::*;
@@ -468,3 +469,32 @@ where
         Ok(self.factorizec_into(UPLO::Upper)?.ln_detc_into())
     }
 }
+
+pub struct PsdCholeskyFactorization {
+    pub pivot: Vec<usize>,
+    pub positive_definite_factorization: CholeskyFactorized<OwnedRepr<f64>>,
+}
+
+pub fn factorize_cholesky_semi_definite(
+    input_mat: &mut ArrayBase<OwnedRepr<f64>, Ix2>,
+) -> Result<PsdCholeskyFactorization> {
+    let uplo = UPLO::Lower;
+    let mut rank: i32 = 0;
+    let pivot: Vec<usize> = cholesky_semi(
+        input_mat.square_layout()?,
+        uplo,
+        input_mat.as_allocated_mut()?,
+        &mut rank,
+    )?
+    .iter()
+    .map(|x| *x as usize)
+    .collect();
+    let cholesky_factors = input_mat.into_triangular(uplo).slice(s![0..rank, 0..rank]);
+    Ok(PsdCholeskyFactorization {
+        pivot: pivot,
+        positive_definite_factorization: CholeskyFactorized::<OwnedRepr<f64>> {
+            factor: replicate(&cholesky_factors),
+            uplo: uplo,
+        },
+    })
+}

src/lapack/cholesky_semi.rs

@@ -0,0 +1,25 @@
+use super::*;
+use crate::{error::*, layout::*};
+
+pub fn cholesky_semi(l: MatrixLayout, uplo: UPLO, a: &mut [f64], rank: &mut i32) -> Result<Vec<i32>> {
+    let (n, _) = l.size();
+    let mut ipiv = vec![0; n as usize];
+    let tol = 1.0e-12_f64;
+    let info = unsafe {
+        lapacke::dpstrf(
+            l.lapacke_layout(),
+            uplo as u8,
+            n,
+            a,
+            l.lda(),
+            &mut ipiv,
+            rank,
+            tol,
+        )
+    };
+    if info < 0 {
+        Err(LinalgError::Lapack { return_code: info })
+    } else {
+        Ok(ipiv)
+    }
+}

src/lapack/mod.rs

@@ -1,6 +1,7 @@
 //! Define traits wrapping LAPACK routines
 
 pub mod cholesky;
+pub mod cholesky_semi;
 pub mod eig;
 pub mod eigh;
 pub mod least_squares;
@@ -14,6 +15,7 @@ pub mod triangular;
 pub mod tridiagonal;
 
 pub use self::cholesky::*;
+pub use self::cholesky_semi::*;
 pub use self::eig::*;
 pub use self::eigh::*;
 pub use self::least_squares::*;