diff --git a/src/SofaCaribou/Python/CMakeLists.txt b/src/SofaCaribou/Python/CMakeLists.txt
index 112c9ebb..a63f95e5 100644
--- a/src/SofaCaribou/Python/CMakeLists.txt
+++ b/src/SofaCaribou/Python/CMakeLists.txt
@@ -41,7 +41,15 @@ set(PYTHON_TEST_FILES
 )
 
 find_package(SofaPython3 REQUIRED)
-find_package(MKL REQUIRED)
+
+if (CARIBOU_WITH_MKL)
+    set(MKL_STATIC ON)
+    if (NOT CARIBOU_WITH_OPENMP)
+        set(MKL_THREADING_VENDOR SEQUENTIAL)
+    endif()
+    find_package(MKL REQUIRED QUIET)
+endif()
+
 caribou_add_python_module(SofaCaribou
                           TARGET_NAME                       ${PROJECT_NAME}
                           TARGET_ALIAS                      Caribou::Python.SofaCaribou
diff --git a/src/SofaCaribou/Python/Solver/ConjugateGradientSolver.h b/src/SofaCaribou/Python/Solver/ConjugateGradientSolver.h
index c02ce1a7..5b0b21f2 100644
--- a/src/SofaCaribou/Python/Solver/ConjugateGradientSolver.h
+++ b/src/SofaCaribou/Python/Solver/ConjugateGradientSolver.h
@@ -13,17 +13,17 @@ namespace SofaCaribou::solver::python {
 template <typename EigenMatrix>
 void bind_ConjugateGradientSolver(pybind11::module & m) {
     namespace py = pybind11;
-    using SOLVER = ConjugateGradientSolver<EigenMatrix>;
-    py::class_<SOLVER, sofa::core::objectmodel::BaseObject, sofapython3::py_shared_ptr<SOLVER>> c(m, "ConjugateGradientSolver");
+    using SolverType = ConjugateGradientSolver<EigenMatrix>;
+    py::class_<SolverType, sofa::core::objectmodel::BaseObject, sofapython3::py_shared_ptr<SolverType>> c(m, "ConjugateGradientSolver");
 
-    c.def("A", [](SOLVER & solver){return solver.A()->matrix();});
+    c.def("A", [](const SolverType & solver){return solver.A()->matrix();});
 
-    c.def("x", [](SOLVER & solver){return solver.x()->vector();});
+    c.def("x", [](const SolverType & solver){return solver.x()->vector();});
 
-    c.def("b", [](SOLVER & solver){return solver.b()->vector();});
+    c.def("b", [](const SolverType & solver){return solver.b()->vector();});
 
 
-    c.def("assemble", [](SOLVER & solver, double m, double b, double k) {
+    c.def("assemble", [](SolverType & solver, double m, double b, double k) {
         sofa::core::MechanicalParams mparams;
         mparams.setMFactor(m);
         mparams.setBFactor(b);
@@ -31,8 +31,8 @@ void bind_ConjugateGradientSolver(pybind11::module & m) {
         solver.assemble(&mparams);
     }, py::arg("m") = static_cast<double>(1), py::arg("b") = static_cast<double>(1), py::arg("k") = static_cast<double>(1));
 
-    sofapython3::PythonFactory::registerType<SOLVER>([](sofa::core::objectmodel::Base* o) {
-        return py::cast(dynamic_cast<SOLVER*>(o));
+    sofapython3::PythonFactory::registerType<SolverType>([](sofa::core::objectmodel::Base* o) {
+        return py::cast(dynamic_cast<SolverType*>(o));
     });
 }
 
diff --git a/src/SofaCaribou/Python/Solver/LDLTSolver.cpp b/src/SofaCaribou/Python/Solver/LDLTSolver.cpp
index 9b774e62..b1ba7ed6 100644
--- a/src/SofaCaribou/Python/Solver/LDLTSolver.cpp
+++ b/src/SofaCaribou/Python/Solver/LDLTSolver.cpp
@@ -7,7 +7,6 @@ namespace SofaCaribou::solver::python {
 
     void addLDLTSolver(py::module & m) {
         bind_LDLTSolver<Eigen::SimplicialLDLT<Eigen::SparseMatrix<FLOATING_POINT_TYPE, Eigen::ColMajor, int>, Eigen::Lower, Eigen::AMDOrdering<int>>>(m);
-
     }
 
 } // namespace SofaCaribou::solver::python
\ No newline at end of file
diff --git a/src/SofaCaribou/Python/Solver/LDLTSolver.h b/src/SofaCaribou/Python/Solver/LDLTSolver.h
index 838ca925..52ff7a6b 100644
--- a/src/SofaCaribou/Python/Solver/LDLTSolver.h
+++ b/src/SofaCaribou/Python/Solver/LDLTSolver.h
@@ -14,17 +14,17 @@ namespace SofaCaribou::solver::python {
     template <typename EigenSolver>
     void bind_LDLTSolver(pybind11::module & m) {
         namespace py = pybind11;
-        using SOLVER = SofaCaribou::solver::LDLTSolver<EigenSolver>;
-        py::class_<SOLVER, sofa::core::objectmodel::BaseObject, sofapython3::py_shared_ptr<SOLVER>> c(m, "LDLTSolver");
+        using SolverType = SofaCaribou::solver::LDLTSolver<EigenSolver>;
+        py::class_<SolverType, sofa::core::objectmodel::BaseObject, sofapython3::py_shared_ptr<SolverType>> c(m, "LDLTSolver");
 
-        c.def("A", [](SOLVER & solver){return solver.A()->matrix();});
+        c.def("A", [](const SolverType & solver){return solver.A()->matrix();});
 
-        c.def("x", [](SOLVER & solver){return solver.x()->vector();});
+        c.def("x", [](const SolverType & solver){return solver.x()->vector();});
 
-        c.def("b", [](SOLVER & solver){return solver.b()->vector();});
+        c.def("b", [](const SolverType & solver){return solver.b()->vector();});
 
 
-        c.def("assemble", [](SOLVER & solver, double m, double b, double k) {
+        c.def("assemble", [](SolverType & solver, double m, double b, double k) {
             sofa::core::MechanicalParams mparams;
             mparams.setMFactor(m);
             mparams.setBFactor(b);
@@ -32,8 +32,8 @@ namespace SofaCaribou::solver::python {
             solver.assemble(&mparams);
         }, py::arg("m") = static_cast<double>(1), py::arg("b") = static_cast<double>(1), py::arg("k") = static_cast<double>(1));
 
-        sofapython3::PythonFactory::registerType<SOLVER>([](sofa::core::objectmodel::Base* o) {
-            return py::cast(dynamic_cast<SOLVER*>(o));
+        sofapython3::PythonFactory::registerType<SolverType>([](sofa::core::objectmodel::Base* o) {
+            return py::cast(dynamic_cast<SolverType*>(o));
         });
     }
 
diff --git a/src/SofaCaribou/Python/Solver/LLTSolver.cpp b/src/SofaCaribou/Python/Solver/LLTSolver.cpp
index fde85921..567d2f74 100644
--- a/src/SofaCaribou/Python/Solver/LLTSolver.cpp
+++ b/src/SofaCaribou/Python/Solver/LLTSolver.cpp
@@ -6,7 +6,6 @@ namespace SofaCaribou::solver::python {
 
     void addLLTSolver(py::module & m) {
         bind_LLTSolver<Eigen::SimplicialLLT<Eigen::SparseMatrix<FLOATING_POINT_TYPE, Eigen::ColMajor, int>, Eigen::Lower, Eigen::AMDOrdering<int>>>(m);
-
     }
 
 } // namespace SofaCaribou::solver::python
\ No newline at end of file
diff --git a/src/SofaCaribou/Python/Solver/LLTSolver.h b/src/SofaCaribou/Python/Solver/LLTSolver.h
index 8f1587ee..263e574d 100644
--- a/src/SofaCaribou/Python/Solver/LLTSolver.h
+++ b/src/SofaCaribou/Python/Solver/LLTSolver.h
@@ -14,16 +14,16 @@ namespace SofaCaribou::solver::python {
     template <typename EigenSolver>
     void bind_LLTSolver(pybind11::module & m) {
         namespace py = pybind11;
-        using SOLVER = LLTSolver<EigenSolver>;
-        py::class_<SOLVER, sofa::core::objectmodel::BaseObject, sofapython3::py_shared_ptr<SOLVER>> c(m, "LLTSolver");
+        using SolverType = LLTSolver<EigenSolver>;
+        py::class_<SolverType, sofa::core::objectmodel::BaseObject, sofapython3::py_shared_ptr<SolverType>> c(m, "LLTSolver");
 
-        c.def("A", [](SOLVER & solver){return solver.A()->matrix();});
+        c.def("A", [](const SolverType & solver){return solver.A()->matrix();});
 
-        c.def("x", [](SOLVER & solver){return solver.x()->vector();});
+        c.def("x", [](const SolverType & solver){return solver.x()->vector();});
 
-        c.def("b", [](SOLVER & solver){return solver.b()->vector();});
+        c.def("b", [](const SolverType & solver){return solver.b()->vector();});
 
-        c.def("assemble", [](SOLVER & solver, double m, double b, double k) {
+        c.def("assemble", [](SolverType & solver, double m, double b, double k) {
             sofa::core::MechanicalParams mparams;
             mparams.setMFactor(m);
             mparams.setBFactor(b);
@@ -31,8 +31,8 @@ namespace SofaCaribou::solver::python {
             solver.assemble(&mparams);
         }, py::arg("m") = static_cast<double>(1), py::arg("b") = static_cast<double>(1), py::arg("k") = static_cast<double>(1));
 
-        sofapython3::PythonFactory::registerType<SOLVER>([](sofa::core::objectmodel::Base* o) {
-            return py::cast(dynamic_cast<SOLVER*>(o));
+        sofapython3::PythonFactory::registerType<SolverType>([](sofa::core::objectmodel::Base* o) {
+            return py::cast(dynamic_cast<SolverType*>(o));
         });
     }
 
diff --git a/src/SofaCaribou/Python/Solver/LUSolver.cpp b/src/SofaCaribou/Python/Solver/LUSolver.cpp
index 0f6241bc..775f376a 100644
--- a/src/SofaCaribou/Python/Solver/LUSolver.cpp
+++ b/src/SofaCaribou/Python/Solver/LUSolver.cpp
@@ -6,7 +6,6 @@ namespace SofaCaribou::solver::python {
 
     void addLUSolver(py::module & m) {
         bind_LUSolver<Eigen::SparseLU<Eigen::SparseMatrix<FLOATING_POINT_TYPE, Eigen::ColMajor, int>, Eigen::AMDOrdering<int>>>(m);
-
     }
 
 } // namespace SofaCaribou::solver::python
\ No newline at end of file
diff --git a/src/SofaCaribou/Python/Solver/LUSolver.h b/src/SofaCaribou/Python/Solver/LUSolver.h
index 845850f1..89cd6580 100644
--- a/src/SofaCaribou/Python/Solver/LUSolver.h
+++ b/src/SofaCaribou/Python/Solver/LUSolver.h
@@ -14,16 +14,16 @@ namespace SofaCaribou::solver::python {
     template <typename EigenSolver>
     void bind_LUSolver(pybind11::module & m) {
         namespace py = pybind11;
-        using SOLVER = LUSolver<EigenSolver>;
-        py::class_<SOLVER, sofa::core::objectmodel::BaseObject, sofapython3::py_shared_ptr<SOLVER>> c(m, "LUSolver");
+        using SolverType = LUSolver<EigenSolver>;
+        py::class_<SolverType, sofa::core::objectmodel::BaseObject, sofapython3::py_shared_ptr<SolverType>> c(m, "LUSolver");
 
-        c.def("A", [](SOLVER & solver){return solver.A()->matrix();});
+        c.def("A", [](const SolverType & solver){return solver.A()->matrix();});
 
-        c.def("x", [](SOLVER & solver){return solver.x()->vector();});
+        c.def("x", [](const SolverType & solver){return solver.x()->vector();});
 
-        c.def("b", [](SOLVER & solver){return solver.b()->vector();});
+        c.def("b", [](const SolverType & solver){return solver.b()->vector();});
 
-        c.def("assemble", [](SOLVER & solver, double m, double b, double k) {
+        c.def("assemble", [](SolverType & solver, double m, double b, double k) {
             sofa::core::MechanicalParams mparams;
             mparams.setMFactor(m);
             mparams.setBFactor(b);
@@ -31,8 +31,8 @@ namespace SofaCaribou::solver::python {
             solver.assemble(&mparams);
         }, py::arg("m") = static_cast<double>(1), py::arg("b") = static_cast<double>(1), py::arg("k") = static_cast<double>(1));
 
-        sofapython3::PythonFactory::registerType<SOLVER>([](sofa::core::objectmodel::Base* o) {
-            return py::cast(dynamic_cast<SOLVER*>(o));
+        sofapython3::PythonFactory::registerType<SolverType>([](sofa::core::objectmodel::Base* o) {
+            return py::cast(dynamic_cast<SolverType*>(o));
         });
     }
 
diff --git a/src/SofaCaribou/SofaCaribouConfig.cmake.in b/src/SofaCaribou/SofaCaribouConfig.cmake.in
index 2a4f2b0a..6765c5fd 100644
--- a/src/SofaCaribou/SofaCaribouConfig.cmake.in
+++ b/src/SofaCaribou/SofaCaribouConfig.cmake.in
@@ -3,6 +3,7 @@
 # OPTIONS
 set(CARIBOU_WITH_OPENMP "@CARIBOU_WITH_OPENMP@")
 set(CARIBOU_WITH_SP3    "@CARIBOU_WITH_SP3@")
+set(CARIBOU_WITH_MKL    "@CARIBOU_WITH_MKL@")
 
 # Caribou
 find_package(Caribou COMPONENTS Algebra Geometry Topology Mechanics REQUIRED)
@@ -13,6 +14,16 @@ if (CARIBOU_WITH_SP3)
 find_package(Caribou COMPONENTS Python REQUIRED)
 endif()
 
+# MKL support
+if (CARIBOU_WITH_MKL)
+    set(MKL_STATIC ON)
+    if (NOT CARIBOU_WITH_OPENMP)
+        set(MKL_THREADING_VENDOR SEQUENTIAL)
+    endif()
+    find_package(MKL REQUIRED QUIET)
+endif()
+
+
 # Sofa's packages
 find_package(SOFA COMPONENTS SofaFramework SofaBaseLinearSolver SofaBaseTopology SofaEigen2Solver QUIET MODULE REQUIRED)