update to faer v0.21

Cargo.toml

@@ -60,7 +60,7 @@ julia = ["sdp", "dep:libc", "dep:num-derive",  "serde", "faer-sparse", "buildinf
 python = ["sdp", "dep:libc", "dep:pyo3", "dep:num-derive", "serde", "faer-sparse", "buildinfo"]
 
 #compile with faer supernodal solver option
-faer-sparse = ["dep:faer", "dep:faer-entity"]
+faer-sparse = ["dep:faer", "dep:faer-traits"]
 
 # -------------------------------
 # SDP configuration
@@ -96,11 +96,11 @@ optional = true
 # -------------------------------
 
 [dependencies.faer]
-version = "0.20"
+version = "0.21"
 optional = true 
 
-[dependencies.faer-entity]
-version = "0.20"
+[dependencies.faer-traits]
+version = "0.21"
 optional = true 
 
 

build.rs

@@ -54,6 +54,7 @@ fn config_python_blas() {
         feature = "sdp-mkl",
         feature = "sdp-r"
     )) {
+        println!("cargo:rustc-cfg=sdp_pyblas");
         printinfo!("Python: compiling with python blas from scipy");
     } else {
         printinfo!("Python: compiling with local blas/lapack libraries");

src/algebra/dense/blas/traits.rs

@@ -6,7 +6,7 @@ cfg_if::cfg_if! {
         // imports via scipy 
         use crate::python::pyblas::*;
     }
-        else {   
+    else {
         // standard imports via blas-lapack-rs crates 
         extern crate blas_src;
         extern crate lapack_src;

src/algebra/floats.rs

@@ -65,8 +65,8 @@ cfg_if::cfg_if! {
 
 cfg_if::cfg_if! {
     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> {}
+        pub trait MaybeFaerFloatT : faer_traits::RealField {}
+        impl<T> MaybeFaerFloatT for T where T: faer_traits::RealField {}
     }
     else {
         pub trait MaybeFaerFloatT {}

src/julia/ClarabelRs/src/ClarabelRs.jl

@@ -22,7 +22,7 @@ module ClarabelRs
         global libpath
         global librust
 
-        result = true
+        result = false
         while !result
             result = Libdl.dlclose(librust)
         end 

src/julia/ClarabelRs/test/test_problem_depot.jl

@@ -1,20 +1,18 @@
 using Convex, Clarabel, ClarabelRs, Test
 
-if false 
+#=
 #tests everything 
-    @testset "ClarabelRs" begin
-    Convex.ProblemDepot.run_tests([r""]; exclude=[r"mip"]) do p
-        solve!(p, Convex.MOI.OptimizerWithAttributes(
-            ClarabelRs.Optimizer, "verbose" => true, 
-            "chordal_decomposition_compact" => true,
-            "chordal_decomposition_enable" => true,
-            "chordal_decomposition_merge_method" => :clique_graph,
-            "chordal_decomposition_complete_dual" => true,
-            ))
-    end
+@testset "ClarabelRs" begin
+Convex.ProblemDepot.run_tests([r""]; exclude=[r"mip"]) do p
+    solve!(p, Convex.MOI.OptimizerWithAttributes(
+        ClarabelRs.Optimizer, "verbose" => true, 
+        "chordal_decomposition_compact" => true,
+        "chordal_decomposition_enable" => true,
+        "chordal_decomposition_merge_method" => :clique_graph,
+        "chordal_decomposition_complete_dual" => true,
+        ))
 end
-
-end 
+=#
 
 @testset "ClarabelRs" begin
     Convex.ProblemDepot.run_tests([r"sdp_quantum_relative"]; exclude=[r"mip",r"sdp_quantum_relative_entropy3_lowrank"]) do p

src/python/impl_default_py.rs

@@ -210,7 +210,9 @@ pub struct PyDefaultSettings {
     #[pyo3(get, set)]
     pub min_terminate_step_length: f64,
 
-    // KKT settings incomplete
+    // KKT settings
+    #[pyo3(get, set)]
+    pub max_threads: u32,
     #[pyo3(get, set)]
     pub direct_kkt_solver: bool,
     #[pyo3(get, set)]
@@ -312,6 +314,7 @@ impl PyDefaultSettings {
             linesearch_backtrack_step: set.linesearch_backtrack_step,
             min_switch_step_length: set.min_switch_step_length,
             min_terminate_step_length: set.min_terminate_step_length,
+            max_threads: set.max_threads,
             direct_kkt_solver: set.direct_kkt_solver,
             direct_solve_method: set.direct_solve_method.clone(),
             static_regularization_enable: set.static_regularization_enable,
@@ -360,6 +363,7 @@ impl PyDefaultSettings {
             linesearch_backtrack_step: self.linesearch_backtrack_step,
             min_switch_step_length: self.min_switch_step_length,
             min_terminate_step_length: self.min_terminate_step_length,
+            max_threads: self.max_threads,
             direct_kkt_solver: self.direct_kkt_solver,
             direct_solve_method: self.direct_solve_method.clone(),
             static_regularization_enable: self.static_regularization_enable,

src/python/io.rs

@@ -76,3 +76,7 @@ pub(crate) fn stdout() -> PythonStdout {
 pub(crate) fn stderr() -> PythonStderr {
     PythonStderr::new()
 }
+
+#[allow(unused_imports)]
+pub(crate) use PythonStderr as Stderr;
+pub(crate) use PythonStdout as Stdout;

src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/faer_ldl.rs

@@ -1,13 +1,29 @@
 #![allow(non_snake_case)]
 
+use faer::dyn_stack::{MemBuffer, MemStack, StackReq};
+// use faer::{
+//     dyn_stack::{GlobalPodBuffer, PodStack, StackReq},
+//     linalg::cholesky::ldlt::compute::LdltRegularization,
+//     sparse::{
+//         linalg::{amd::Control, cholesky::*, SupernodalThreshold},
+//         SparseColMatRef, SymbolicSparseColMatRef,
+//     },
+//     Conj, Parallelism, Side,
+// };
 use faer::{
-    dyn_stack::{GlobalPodBuffer, PodStack, StackReq},
-    linalg::cholesky::ldlt_diagonal::compute::LdltRegularization,
+    linalg::cholesky::ldlt::factor::{LdltParams, LdltRegularization},
     sparse::{
-        linalg::{amd::Control, cholesky::*, SupernodalThreshold},
+        linalg::{
+            amd::Control,
+            cholesky::{
+                factorize_symbolic_cholesky, CholeskySymbolicParams, LdltRef, SymbolicCholesky,
+                SymmetricOrdering,
+            },
+            SupernodalThreshold,
+        },
         SparseColMatRef, SymbolicSparseColMatRef,
     },
-    Conj, Parallelism, Side,
+    Conj, MatMut, Par, Side, Spec,
 };
 
 use crate::algebra::*;
@@ -56,29 +72,39 @@ pub struct FaerDirectLDLSolver<T: FloatT> {
 
     // regularizer parameters captured from settings
     regularizer_params: FaerLDLRegularizerParams<T>,
+    ldlt_params: Spec<LdltParams, T>,
 
     // symbolic + numeric cholesky data
     symbolic_cholesky: SymbolicCholesky<usize>,
     ld_vals: Vec<T>,
 
     // workspace for faer factor/solve calls
-    work: GlobalPodBuffer,
+    work: MemBuffer,
 
-    parallelism: Parallelism<'static>,
+    parallelism: Par,
 }
 
 impl<T> FaerDirectLDLSolver<T>
 where
     T: FloatT,
 {
+    pub fn nthreads_from_settings(setting: usize) -> usize {
+        faer::utils::thread::parallelism_degree(faer::Par::rayon(setting))
+    }
+
     pub fn new(KKT: &CscMatrix<T>, Dsigns: &[i8], settings: &CoreSettings<T>) -> Self {
         assert!(KKT.is_square(), "KKT matrix is not square");
 
         // -----------------------------
 
-        // Rayon(0) here is equivalent to rayon::current_num_threads()
-        // PJG: this should be a user-settable option
-        let parallelism = Parallelism::Rayon(0);
+        // Par::rayon(0) here is equivalent to rayon::current_num_threads()
+        let parallelism = {
+            match settings.max_threads {
+                0 => Par::rayon(0),
+                1 => Par::Seq,
+                _ => Par::rayon(settings.max_threads as usize),
+            }
+        };
 
         // manually compute an AMD ordering for the KKT matrix
         // and permute it to match the ordering used in QDLDL
@@ -125,21 +151,22 @@ where
         )
         .unwrap();
 
-        let ld_vals = vec![T::zero(); symbolic_cholesky.len_values()];
+        let ld_vals = vec![T::zero(); symbolic_cholesky.len_val()];
 
         let regularizer_params = FaerLDLRegularizerParams {
             regularize: settings.dynamic_regularization_enable,
             delta: settings.dynamic_regularization_delta,
             eps: settings.dynamic_regularization_eps,
         };
 
+        let ldlt_params: Spec<LdltParams, T> = Default::default();
+
         // Required workspace for faer factor and solve
-        let req_factor = symbolic_cholesky
-            .factorize_numeric_ldlt_req::<f64>(true, parallelism)
-            .unwrap();
-        let req_solve = symbolic_cholesky.solve_in_place_req::<f64>(1).unwrap();
-        let req = StackReq::any_of([req_factor, req_solve]);
-        let work = GlobalPodBuffer::new(req);
+        let req_factor =
+            symbolic_cholesky.factorize_numeric_ldlt_scratch::<T>(parallelism, ldlt_params);
+        let req_solve = symbolic_cholesky.solve_in_place_scratch::<T>(1, parallelism); // 1 is the number of RHS
+        let req = StackReq::any_of(&[req_factor, req_solve]);
+        let work = MemBuffer::new(req);
 
         let bperm = vec![T::zero(); perm_kkt.n];
 
@@ -151,6 +178,7 @@ where
             perm_dsigns,
             bperm,
             regularizer_params,
+            ldlt_params,
             symbolic_cholesky,
             ld_vals,
             work,
@@ -161,7 +189,7 @@ where
 
 impl<T> DirectLDLSolver<T> for FaerDirectLDLSolver<T>
 where
-    T: FloatT + faer::Entity + faer::SimpleEntity,
+    T: FloatT,
 {
     fn update_values(&mut self, index: &[usize], values: &[T]) {
         // PJG: this is replicating the update_values function in qdldl
@@ -198,18 +226,14 @@ where
         // NB: faer solves in place.  Permute b to match the ordering used internally
         permute(&mut self.bperm, b, &self.perm);
 
-        let rhs = faer::mat::from_column_major_slice_mut::<T, usize, usize>(
-            &mut self.bperm[0..],
-            b.len(),
-            1,
-        );
+        let rhs = MatMut::from_column_major_slice_mut(&mut self.bperm[0..], b.len(), 1);
         let ldlt = LdltRef::new(&self.symbolic_cholesky, self.ld_vals.as_slice());
 
         ldlt.solve_in_place_with_conj(
             Conj::No,
             rhs,
             self.parallelism,
-            PodStack::new(&mut self.work),
+            MemStack::new(&mut self.work),
         );
 
         // bperm is now the solution, permute it back to the original ordering
@@ -230,16 +254,17 @@ where
 
         let regularizer = self.regularizer_params.to_faer(&self.perm_dsigns);
 
-        self.symbolic_cholesky.factorize_numeric_ldlt(
-            self.ld_vals.as_mut_slice(),
-            a,
-            Side::Upper,
-            regularizer,
-            self.parallelism,
-            PodStack::new(&mut self.work),
-        );
-
-        true // assume success for experimental purposes
+        self.symbolic_cholesky
+            .factorize_numeric_ldlt(
+                self.ld_vals.as_mut_slice(),
+                a,
+                Side::Upper,
+                regularizer,
+                self.parallelism,
+                MemStack::new(&mut self.work),
+                self.ldlt_params,
+            )
+            .is_ok() // PJG: convert to bool for consistency with qdldl.   Should really return Result here and elsewhere
     }
 
     fn required_matrix_shape() -> MatrixTriangle {

src/solver/implementations/default/info_print.rs

@@ -172,12 +172,14 @@ fn _print_settings<T: FloatT>(settings: &DefaultSettings<T>) -> std::io::Result<
     writeln!(out, "settings:")?;
 
     if set.direct_kkt_solver {
-        writeln!(
+        write!(
             out,
             "  linear algebra: direct / {}, precision: {} bit",
             set.direct_solve_method,
             _get_precision_string::<T>()
         )?;
+        print_nthreads(&mut out, settings)?;
+        writeln!(out)?;
     }
 
     let time_lim_str = {
@@ -246,6 +248,28 @@ fn _print_settings<T: FloatT>(settings: &DefaultSettings<T>) -> std::io::Result<
     std::io::Result::Ok(())
 }
 
+#[allow(unused_variables)] //out is unused if faer-sparse is not enabled
+fn print_nthreads<T: FloatT>(
+    out: &mut stdio::Stdout,
+    settings: &DefaultSettings<T>,
+) -> std::io::Result<()> {
+    match settings.direct_solve_method.as_str() {
+        #[cfg(feature = "faer-sparse")]
+        "faer" => {
+            let nthreads =
+                crate::solver::core::kktsolvers::direct::ldlsolvers::faer_ldl::FaerDirectLDLSolver::<
+                    T,
+                >::nthreads_from_settings(settings.max_threads as usize);
+            if nthreads == 1 {
+                write!(out, " ({nthreads} thread) ")
+            } else {
+                write!(out, " ({nthreads} threads) ")
+            }
+        }
+        _ => std::io::Result::Ok(()),
+    }
+}
+
 #[cfg(feature = "sdp")]
 fn print_chordal_decomposition<T: FloatT>(
     chordal_info: &ChordalInfo<T>,

src/solver/implementations/default/settings.rs

@@ -108,6 +108,11 @@ pub struct DefaultSettings<T: FloatT> {
     #[builder(default = "(1e-4).as_T()")]
     pub min_terminate_step_length: T,
 
+    ///maximum solver threads for multithreaded KKT solvers
+    ///choosing 0 lets the solver choose for itself
+    #[builder(default = "0")]
+    pub max_threads: u32,
+
     ///use a direct linear solver method (required true)
     #[builder(default = "true")]
     pub direct_kkt_solver: bool,

src/stdio/mod.rs

@@ -1,10 +1,10 @@
-#![allow(unused_imports)]
+#[cfg(not(feature = "python"))]
+#[allow(unused_imports)]
+pub(crate) use std::io::{stderr, stdout, Stdout};
 
 // configure python specific stdout and stdin streams
 // when compiled with the python feature.   This avoids
 // problems when running within python notebooks etc.
 #[cfg(feature = "python")]
-pub(crate) use crate::python::io::{stderr, stdout};
-
-#[cfg(not(feature = "python"))]
-pub(crate) use std::io::{stderr, stdout};
+#[allow(unused_imports)]
+pub(crate) use crate::python::io::{stderr, stdout, Stdout};