supernodal LDL wrapper (#112)

.github/workflows/ci.yml

@@ -28,7 +28,7 @@ jobs:
 
     - name: Generate code coverage 
       run: |
-        cargo tarpaulin --out xml --features sdp-accelerate,serde --exclude-files "src/python/*,src/julia/*"
+        cargo tarpaulin --out xml --features sdp-accelerate,serde,faer-sparse --exclude-files "src/python/*,src/julia/*"
 
     - name: Upload to codecov.io
       uses: codecov/codecov-action@v4

Cargo.toml

@@ -56,6 +56,8 @@ julia = ["sdp", "dep:libc", "dep:num-derive",  "serde"]
 python = ["sdp", "dep:libc", "dep:pyo3", "dep:num-derive", "serde"]
 
 
+#compile with faer supernodal solver option
+faer-sparse = ["dep:faer", "dep:faer-entity"]
 
 # -------------------------------
 # SDP configuration
@@ -81,6 +83,14 @@ optional = true
 version = "2.2"
 optional = true
 
+[dependencies.faer]
+version = "0.19"
+optional = true 
+
+[dependencies.faer-entity]
+version = "0.19"
+optional = true 
+
 
 # -------------------------------
 # examples
@@ -115,13 +125,19 @@ name = "sdp"
 path = "examples/rust/example_sdp.rs"
 required-features = ["sdp"]
 
+[[example]]
+name = "box_faer"
+path = "examples/rust/example_box_faer.rs"
+required-features = ["faer-sparse"]
+
 [[example]]
 name = "json"
 path = "examples/rust/example_json.rs"
 required-features = ["serde"]
 
 
 
+
 # -------------------------------
 # custom build profiles 
 # -------------------------------

examples/rust/example_box.rs

@@ -3,7 +3,7 @@ use clarabel::algebra::*;
 use clarabel::solver::*;
 
 fn problem_data() -> (CscMatrix<f64>, Vec<f64>, CscMatrix<f64>, Vec<f64>) {
-    let n = 20000;
+    let n = 2000000;
 
     let P = CscMatrix::identity(n);
 

examples/rust/example_box_faer.rs

@@ -0,0 +1,39 @@
+#![allow(non_snake_case)]
+use clarabel::algebra::*;
+use clarabel::solver::*;
+
+fn problem_data() -> (CscMatrix<f64>, Vec<f64>, CscMatrix<f64>, Vec<f64>) {
+    let n = 2000000;
+
+    let P = CscMatrix::identity(n);
+
+    // construct A = [I; -I]
+    let I1 = CscMatrix::<f64>::identity(n);
+    let mut I2 = CscMatrix::<f64>::identity(n);
+    I2.negate();
+
+    let A = CscMatrix::vcat(&I1, &I2);
+
+    let q = vec![1.; n];
+    let b = vec![1.; 2 * n];
+
+    (P, q, A, b)
+}
+
+fn main() {
+    let (P, q, A, b) = problem_data();
+
+    let cones = [NonnegativeConeT(b.len())];
+
+    let settings = DefaultSettingsBuilder::default()
+        .direct_solve_method("faer".to_owned())
+        .equilibrate_enable(true)
+        .max_iter(50)
+        .verbose(true)
+        .build()
+        .unwrap();
+
+    let mut solver = DefaultSolver::new(&P, &q, &A, &b, &cones, settings);
+
+    solver.solve();
+}

src/algebra/floats.rs

@@ -40,53 +40,50 @@ impl<T> CoreFloatT for T where
 {
 }
 
+// additional traits that add bounds to CoreFloatT
+// when optional dependencies are enabled
+
 // if "sdp" is enabled, we must add an additional trait
 // trait bound to restrict compilation for f32/f64 types
 // since there is no BLAS support otherwise
 
-// Define the documentation string as a macro so that FloatT gets documentation
-// regardless of whether the "sdp" feature is enabled.
-macro_rules! floatT_doc_header {
-    () => {
-        r#"Main trait for floating point types used in the Clarabel solver."#
-    };
-}
-macro_rules! floatT_doc_long {
-    () => {
-        r#"All floating point calculations in Clarabel are represented internally on values 
-         implementing the `FloatT` trait, with implementations provided only for f32 and f64 
-         native types when compiled with BLAS/LAPACK support for SDPs. If SDP support is not 
-         enabled then it should be possible to compile Clarabel to support any any other 
-         floating point type provided that it satisfies the trait bounds of `CoreFloatT`. 
-        \
-        \
-         `FloatT` relies on [`num_traits`](num_traits) for most of its constituent trait bounds."#
-    };
-}
 cfg_if::cfg_if! {
-    if #[cfg(not(feature="sdp"))] {
-    #[doc = floatT_doc_header!()]
-    ///
-    #[doc = floatT_doc_long!()]
-        pub trait FloatT: CoreFloatT {}
-    } else{
-        #[doc = floatT_doc_header!()]
-        ///
-        #[doc = floatT_doc_long!()]
-        ///
-        /// The trait bound `BlasFloatT` is only enforced when compiling with the `sdp` feature.
-        pub trait FloatT: CoreFloatT + BlasFloatT {}
+    if #[cfg(feature="sdp")] {
+        pub trait MaybeBlasFloatT : BlasFloatT {}
+        impl<T> MaybeBlasFloatT for T where T: BlasFloatT {}
+    }
+    else {
+        pub trait MaybeBlasFloatT {}
+        impl<T> MaybeBlasFloatT for T {}
     }
 }
 
+// if "faer" is enabled, we must add an additional bound
+// to restrict compilation to Reals implementing SimpleEntity
+
 cfg_if::cfg_if! {
-    if #[cfg(feature="sdp")] {
-        impl<T> FloatT for T where T: CoreFloatT + BlasFloatT {}
-    } else{
-        impl<T> FloatT for T where T: CoreFloatT {}
+    if #[cfg(feature="faer-sparse")] {
+        pub trait MaybeFaerFloatT : faer::SimpleEntity + faer::ComplexField<Real=Self> {}
+        impl<T> MaybeFaerFloatT for T where T: faer::SimpleEntity + faer::ComplexField<Real=Self> {}
+    }
+    else {
+        pub trait MaybeFaerFloatT {}
+        impl<T> MaybeFaerFloatT for T {}
     }
 }
 
+/// Main trait for floating point types used in the Clarabel solver.
+///
+/// All floating point calculations in Clarabel are represented internally on values
+/// implementing the `FloatT` trait, with implementations provided only for f32 and f64
+/// native types when compiled with BLAS/LAPACK support for SDPs. If SDP support is not
+/// enabled then it should be possible to compile Clarabel to support any any other
+/// floating point type provided that it satisfies the trait bounds of `CoreFloatT`.
+///
+/// `FloatT` relies on [`num_traits`](num_traits) for most of its constituent trait bounds.
+pub trait FloatT: CoreFloatT + MaybeBlasFloatT + MaybeFaerFloatT {}
+impl<T> FloatT for T where T: CoreFloatT + MaybeBlasFloatT + MaybeFaerFloatT {}
+
 /// Trait for convering Rust primitives to [`FloatT`](crate::algebra::FloatT)
 ///
 /// This convenience trait is implemented on f32/64 and u32/64.  This trait

src/algebra/tests/vector.rs

@@ -140,6 +140,9 @@ fn test_norm_scaled() {
 fn test_norm_inf() {
     let x = [-3., 4., -12.];
     assert_eq!(x.norm_inf(), 12.);
+
+    let x = [-3., f64::NAN, -12.];
+    assert!(x.norm_inf().is_nan());
 }
 
 #[test]

src/algebra/utils.rs

@@ -6,9 +6,9 @@
 // which serves as a vectorized version of the std::iter::position
 // returning indices of *all* elements satisfying a predicate
 
+use crate::qdldl;
 use num_traits::Num;
 use std::cmp::Ordering;
-use std::iter::zip;
 
 #[cfg_attr(not(sdp), allow(dead_code))]
 pub(crate) trait PositionAll<T>: Iterator<Item = T> {
@@ -34,12 +34,12 @@ where
 
 // permutation and inverse permutation
 pub(crate) fn permute<T: Copy>(x: &mut [T], b: &[T], p: &[usize]) {
-    zip(p, x).for_each(|(p, x)| *x = b[*p]);
+    qdldl::permute(x, b, p);
 }
 
 #[allow(dead_code)]
 pub(crate) fn ipermute<T: Copy>(x: &mut [T], b: &[T], p: &[usize]) {
-    zip(p, b).for_each(|(p, b)| x[*p] = *b);
+    qdldl::ipermute(x, b, p);
 }
 
 // Construct an inverse permutation from a permutation

src/algebra/vecmath.rs

@@ -117,11 +117,11 @@ impl<T: FloatT> VectorMath<T> for [T] {
     // Returns infinity norm
     fn norm_inf(&self) -> T {
         let mut out = T::zero();
-        for v in self.iter().map(|v| v.abs()) {
+        for &v in self {
             if v.is_nan() {
                 return T::nan();
             }
-            out = if v > out { v } else { out };
+            out = T::max(out, v.abs());
         }
         out
     }

src/qdldl/qdldl.rs

@@ -101,7 +101,7 @@ where
 
         // permute b
         let tmp = &mut self.workspace.fwork;
-        _permute(tmp, b, &self.perm);
+        permute(tmp, b, &self.perm);
 
         //solve in place with tmp as permuted RHS
         _solve(
@@ -113,7 +113,7 @@ where
         );
 
         // inverse permutation to put unpermuted soln in b
-        _ipermute(b, tmp, &self.perm);
+        ipermute(b, tmp, &self.perm);
     }
 
     pub fn update_values(&mut self, indices: &[usize], values: &[T]) {
@@ -195,20 +195,20 @@ fn _qdldl_new<T: FloatT>(
         iperm = _invperm(&_perm)?;
         perm = _perm;
     } else {
-        (perm, iperm) = _get_amd_ordering(Ain, opts.amd_dense_scale);
+        (perm, iperm) = get_amd_ordering(Ain, opts.amd_dense_scale);
     }
 
     //permute to (another) upper triangular matrix and store the
     //index mapping the input's entries to the permutation's entries
-    let (A, AtoPAPt) = _permute_symmetric(Ain, &iperm);
+    let (A, AtoPAPt) = permute_symmetric(Ain, &iperm);
 
     // handle the (possibly permuted) vector of
     // diagonal D signs if one was specified.  Otherwise
     // otherwise all signs are positive
     let mut Dsigns = vec![1_i8; n];
     if let Some(ds) = opts.Dsigns {
         Dsigns = vec![1_i8; n];
-        _permute(&mut Dsigns, &ds, &perm);
+        permute(&mut Dsigns, &ds, &perm);
     }
 
     // allocate workspace
@@ -452,8 +452,8 @@ fn _factor_inner<T: FloatT>(
     Lp[0] = 0;
     let mut acc = 0;
     for (Lp, Lnz) in zip(&mut Lp[1..], Lnz) {
-        *Lp = acc + Lnz;
-        acc = *Lp;
+        acc += Lnz;
+        *Lp = acc;
     }
 
     //  set all y_idx to be 'unused' initially
@@ -681,15 +681,38 @@ fn _ltsolve_unsafe<T: FloatT>(Lp: &[usize], Li: &[usize], Lx: &[T], x: &mut [T])
     }
 }
 
+// Solves D(L+I)'x = b, with x replacing b.  Unchecked version.
+fn _dltsolve_unsafe<T: FloatT>(Lp: &[usize], Li: &[usize], Lx: &[T], Dinv: &[T], x: &mut [T]) {
+    unsafe {
+        for i in (0..x.len()).rev() {
+            let mut s = T::zero();
+            let f = *Lp.get_unchecked(i);
+            let l = *Lp.get_unchecked(i + 1);
+            for (&Lxj, &Lij) in zip(&Lx[f..l], &Li[f..l]) {
+                s += Lxj * (*x.get_unchecked(Lij));
+            }
+
+            let xi = x.get_unchecked_mut(i);
+            *xi *= *Dinv.get_unchecked(i);
+            *xi -= s;
+        }
+    }
+}
+
 // Solves Ax = b where A has given LDL factors, with x replacing b
 fn _solve<T: FloatT>(Lp: &[usize], Li: &[usize], Lx: &[T], Dinv: &[T], b: &mut [T]) {
     // We call the `unsafe`d version of the forward and backward substitution
     // functions here, since the matrix data should be well posed and x of
     // compatible dimensions.   For super safety or debugging purposes, there
     // are also `safe` versions implemented above.
     _lsolve_unsafe(Lp, Li, Lx, b);
-    zip(b.iter_mut(), Dinv).for_each(|(b, d)| *b *= *d);
-    _ltsolve_unsafe(Lp, Li, Lx, b);
+
+    // combined D and L^T solve in one pass
+    _dltsolve_unsafe(Lp, Li, Lx, Dinv, b);
+
+    // in separate passes for reference
+    //zip(b.iter_mut(), Dinv).for_each(|(b, d)| *b *= *d);
+    //_ltsolve_unsafe(Lp, Li, Lx, b);
 }
 
 // Construct an inverse permutation from a permutation
@@ -706,21 +729,32 @@ fn _invperm(p: &[usize]) -> Result<Vec<usize>, QDLDLError> {
     Ok(b)
 }
 
-// internal permutation and inverse permutation
-// functions that require no memory allocations
+// permutation and inverse permutation
+// functions that require no allocation
+// p must be a valid permutation vector
+// in both cases for safety
 
-fn _permute<T: Copy>(x: &mut [T], b: &[T], p: &[usize]) {
-    zip(p, x).for_each(|(p, x)| *x = b[*p]);
+pub(crate) fn permute<T: Copy>(x: &mut [T], b: &[T], p: &[usize]) {
+    debug_assert!(*p.iter().max().unwrap_or(&0) < x.len());
+    unsafe {
+        zip(p, x).for_each(|(p, x)| *x = *b.get_unchecked(*p));
+    }
 }
 
-fn _ipermute<T: Copy>(x: &mut [T], b: &[T], p: &[usize]) {
-    zip(p, b).for_each(|(p, b)| x[*p] = *b);
+pub(crate) fn ipermute<T: Copy>(x: &mut [T], b: &[T], p: &[usize]) {
+    debug_assert!(*p.iter().max().unwrap_or(&0) < x.len());
+    unsafe {
+        zip(p, b).for_each(|(p, b)| *x.get_unchecked_mut(*p) = *b);
+    }
 }
 
 // Given a sparse symmetric matrix `A` (with only upper triangular entries), return
 // permuted sparse symmetric matrix `P` (also only upper triangular) given the
 // inverse permutation vector `iperm`."
-fn _permute_symmetric<T: FloatT>(A: &CscMatrix<T>, iperm: &[usize]) -> (CscMatrix<T>, Vec<usize>) {
+pub(crate) fn permute_symmetric<T: FloatT>(
+    A: &CscMatrix<T>,
+    iperm: &[usize],
+) -> (CscMatrix<T>, Vec<usize>) {
     // perform a number of argument checks
     let (_m, n) = A.size();
     let mut P = CscMatrix::<T>::spalloc((n, n), A.nnz());
@@ -816,7 +850,7 @@ fn _permute_symmetric_inner<T: FloatT>(
     }
 }
 
-fn _get_amd_ordering<T: FloatT>(
+pub(crate) fn get_amd_ordering<T: FloatT>(
     A: &CscMatrix<T>,
     amd_dense_scale: f64,
 ) -> (Vec<usize>, Vec<usize>) {