Add pardiso settings (#61)

goulart-paul · web-flow · commit a16d445e3c96 · 2025-05-05T19:29:51.000+01:00
Add pardiso settings in support of oxfordcontrol/Clarabel.rs#181
diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
@@ -47,16 +47,20 @@ jobs:
       run: |
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_expcone
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_lp
-        valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_callback
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_powcone
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_genpowcone
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_qp
+        valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_qp_f32
+        valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_qp_f64
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_socp
+        valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/c/example_callback
+
 
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/cpp/cpp_example_expcone
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/cpp/cpp_example_lp
-        valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/cpp/cpp_example_callback
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/cpp/cpp_example_powcone
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/cpp/cpp_example_genpowcone
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/cpp/cpp_example_qp
         valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/cpp/cpp_example_socp
+        valgrind --leak-check=full --show-leak-kinds=all --track-origins=yes --verbose --error-exitcode=1 ./examples/cpp/cpp_example_callback
+
diff --git a/Clarabel.rs b/Clarabel.rs
@@ -1 +1 @@
-Subproject commit 3ca4b3b3b26a43f6551ed3f4b0bd5530ef51a8c7
+Subproject commit 6292d0a70f4ad44bbb7262c7dbf9b9518c31f6fa
diff --git a/examples/c/CMakeLists.txt b/examples/c/CMakeLists.txt
@@ -12,6 +12,7 @@ set(C_EXAMPLES
   example_json
   example_print_stream
   example_callback
+  example_pardiso_mkl
 )
 
 # Compile utils.c as a library
diff --git a/examples/c/example_pardiso_mkl.c b/examples/c/example_pardiso_mkl.c
@@ -0,0 +1,95 @@
+// To use float instead of double:
+// #define CLARABEL_USE_FLOAT
+
+#include "utils.h"
+#include <stdio.h>
+#include <Clarabel.h>
+
+int main(void)
+{
+
+    #ifndef FEATURE_PARDISO_MKL
+        printf("This example requires MKL Pardiso support.\n");
+        return 1;
+#else
+
+    // QP Example
+
+    /* From dense matrix:
+     * [[6., 0.],
+     *  [0., 4.]]
+     */
+    ClarabelCscMatrix P;
+    clarabel_CscMatrix_init(
+        &P,
+        2,                          // row
+        2,                          // col
+        (uintptr_t[]){ 0, 1, 2 },   // colptr
+        (uintptr_t[]){ 0, 1 },      // rowval
+        (ClarabelFloat[]){ 6., 4. } // nzval
+    );
+
+    ClarabelFloat q[2] = { -1., -4. };
+
+    /* From dense matrix:
+     * [[ 1., -2.], // <-- LHS of equality constraint (lower bound)
+     *  [ 1.,  0.], // <-- LHS of inequality constraint (upper bound)
+     *  [ 0.,  1.], // <-- LHS of inequality constraint (upper bound)
+     *  [-1.,  0.], // <-- LHS of inequality constraint (lower bound)
+     *  [ 0., -1.]] // <-- LHS of inequality constraint (lower bound)
+     */
+    ClarabelCscMatrix A;
+    clarabel_CscMatrix_init(
+        &A,
+        5,                                             // row
+        2,                                             // col
+        (uintptr_t[]){ 0, 3, 6 },                      // colptr
+        (uintptr_t[]){ 0, 1, 3, 0, 2, 4 },             // rowval
+        (ClarabelFloat[]){ 1., 1., -1., -2., 1., -1. } // nzval
+    );
+
+    ClarabelFloat b[5] = { 0., 1., 1., 1., 1. };
+
+    ClarabelSupportedConeT cones[2] = { ClarabelZeroConeT(1), ClarabelNonnegativeConeT(4) };
+
+    // Settings
+    ClarabelDefaultSettings settings = clarabel_DefaultSettings_default();
+    //settings.direct_solve_method = PARDISO_PANUA;    
+    settings.direct_solve_method = PARDISO_MKL;  
+
+    // turn on verbose output and change some settings 
+    settings.pardiso_verbose = true;
+    settings.pardiso_iparm[1] = 0; // tries minimum degree ordering instead of default 
+
+    // Build solver
+    ClarabelDefaultSolver *solver = clarabel_DefaultSolver_new(
+        &P, // P
+        q,  // q
+        &A, // A
+        b,  // b
+        2,  // n_cones
+        cones, &settings
+    );
+
+    // Solve
+    clarabel_DefaultSolver_solve(solver);
+
+    // Get solution
+    ClarabelDefaultSolution solution = clarabel_DefaultSolver_solution(solver);
+    print_solution(&solution);
+
+    // Get some detailed solve information
+    ClarabelDefaultInfo_f64 info = clarabel_DefaultSolver_info(solver);
+    printf("primal residual = %e\n", info.res_primal);
+    printf("dual residual   = %e\n", info.res_dual);
+    printf("# of threads    = %d\n", info.linsolver.threads);
+    printf("KKT nonzeros    = %d\n", info.linsolver.nnzA);
+    printf("factor nonzeros = %d\n", info.linsolver.nnzL);
+
+    // Free the matrices and the solver
+    clarabel_DefaultSolver_free(solver);
+
+    return 0;
+
+#endif
+}
diff --git a/examples/cpp/CMakeLists.txt b/examples/cpp/CMakeLists.txt
@@ -7,11 +7,11 @@ set(CPP_EXAMPLES
   example_qp
   example_socp
   example_sdp
-  example_faer
-  example_mkl_pardiso
   example_json
   example_print_stream
   example_callback
+  example_faer
+  example_pardiso_mkl
 )
 
 # Define an executable target for each example
diff --git a/examples/cpp/example_pardiso_mkl.cpp b/examples/cpp/example_pardiso_mkl.cpp
@@ -54,7 +54,13 @@ int main(void)
 
     // Settings
     DefaultSettings<double> settings = DefaultSettings<double>::default_settings();
-    settings.direct_solve_method = ClarabelDirectSolveMethods::PARDISO_MKL;
+
+    //settings.direct_solve_method = ClarabelDirectSolveMethods::PARDISO_PANUA;    
+    settings.direct_solve_method = ClarabelDirectSolveMethods::PARDISO_MKL;    
+
+    settings.pardiso_verbose = true;
+    settings.pardiso_iparm[1] = 0; // tries minimum degree ordering instead of default 
+
     //set the number of threads
     settings.max_threads = 16;
 
@@ -68,6 +74,14 @@ int main(void)
     DefaultSolution<double> solution = solver.solution();
     utils::print_solution(solution);
 
+    DefaultInfo<double> info = solver.info();
+
+    printf("primal residual = %e\n", info.res_primal);
+    printf("dual residual   = %e\n", info.res_dual);
+    printf("# of threads    = %d\n", info.linsolver.threads);
+    printf("KKT nonzeros    = %d\n", info.linsolver.nnzA);
+    printf("factor nonzeros = %d\n", info.linsolver.nnzL);
+
     return 0;
 
 #endif // FEATURE_FAER_SPARSE
diff --git a/include/c/DefaultSettings.h b/include/c/DefaultSettings.h
@@ -78,6 +78,10 @@ typedef struct ClarabelDefaultSettings_f64
         bool chordal_decomposition_compact;
         bool chordal_decomposition_complete_dual;
     #endif
+    #ifdef FEATURE_PARDISO_ANY
+        int32_t pardiso_iparm[64];
+        bool pardiso_verbose;
+    #endif 
 } ClarabelDefaultSettings_f64;
 
 typedef struct ClarabelDefaultSettings_f32
@@ -126,6 +130,10 @@ typedef struct ClarabelDefaultSettings_f32
         bool chordal_decomposition_compact;
         bool chordal_decomposition_complete_dual;
     #endif
+    #ifdef FEATURE_PARDISO_ANY
+        int32_t pardiso_iparm[64];
+        bool pardiso_verbose;
+#endif 
 } ClarabelDefaultSettings_f32;
 
 #ifdef CLARABEL_USE_FLOAT
diff --git a/include/cpp/DefaultSettings.h b/include/cpp/DefaultSettings.h
@@ -79,7 +79,10 @@ struct DefaultSettings
         bool chordal_decomposition_compact;
         bool chordal_decomposition_complete_dual;
     #endif
-    
+    #ifdef FEATURE_PARDISO_ANY
+        int32_t pardiso_iparm[64];
+        bool pardiso_verbose;
+    #endif 
 
     static DefaultSettings<T> default_settings();
 };
@@ -355,6 +358,22 @@ class DefaultSettingsBuilder
             return *this;
         }
     #endif
+    #ifdef FEATURE_PARDISO_ANY
+        DefaultSettingsBuilder<T> &pardiso_iparm(int32_t pardiso_iparm[64])
+        {
+            for (int i = 0; i < 64; i++)
+            {
+                settings.pardiso_iparm[i] = pardiso_iparm[i];
+            }
+            return *this;
+        }
+
+        DefaultSettingsBuilder<T> &pardiso_verbose(bool pardiso_verbose)
+        {
+            settings.pardiso_verbose = pardiso_verbose;
+            return *this;
+        }
+    #endif
 
 };
 
diff --git a/rust_wrapper/CMakeLists.txt b/rust_wrapper/CMakeLists.txt
@@ -57,6 +57,7 @@ if(CLARABEL_FEATURE_PARDISO_MKL)
     # Define the FEATURE_PARDISO_MKL flag for all targets that link against clarabel_c
     target_compile_definitions(libclarabel_c_static INTERFACE FEATURE_PARDISO_MKL)
     target_compile_definitions(libclarabel_c_shared INTERFACE FEATURE_PARDISO_MKL)
+    
 endif()
 
 # PARDISO_PANUA feature flag
@@ -70,6 +71,11 @@ if(CLARABEL_FEATURE_PARDISO_PANUA)
     target_compile_definitions(libclarabel_c_shared INTERFACE FEATURE_PARDISO_PANUA)
 endif()
 
+if(CLARABEL_FEATURE_PARDISO_MKL OR CLARABEL_FEATURE_PARDISO_PANUA)
+    target_compile_definitions(libclarabel_c_static INTERFACE FEATURE_PARDISO_ANY)
+    target_compile_definitions(libclarabel_c_shared INTERFACE FEATURE_PARDISO_ANY)
+endif()
+
 
 # SERDE feature flag
 if(CLARABEL_FEATURE_SERDE)

Original file line number	Diff line number	Diff line change
`@@ -12,6 +12,7 @@ set(C_EXAMPLES`
`12`	`12`	`example_json`
`13`	`13`	`example_print_stream`
`14`	`14`	`example_callback`
	`15`	`+ example_pardiso_mkl`
`15`	`16`	`)`
`16`	`17`
`17`	`18`	`# Compile utils.c as a library`