diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 636b860..31568d5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,80 +1,81 @@
 ################## Tests ########################################
 
 include(CatchTest)
 
 
 set(PROJECT_TEST_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 include_directories(${PROJECT_TEST_DIR})
 
 add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND}) # make check is also valid
 
 
 # Catch test
 #===========
 
 #common
 
 set( COMMON_TEST_DIR common )
 
 set(SPECMICP_COMMON_TEST_FILES
     ${COMMON_TEST_DIR}/test_common.cpp
     ${COMMON_TEST_DIR}/units.cpp
     ${COMMON_TEST_DIR}/laws.cpp
     ${COMMON_TEST_DIR}/misc.cpp
     ${COMMON_TEST_DIR}/value_checker.cpp
+    ${COMMON_TEST_DIR}/sparse_solvers.cpp
 )
 
 if (HDF5_FOUND)
     list(APPEND SPECMICP_COMMON_TEST_FILES
         ${COMMON_TEST_DIR}/hdf5_eigen.cpp
     )
 
     include_directories(HDF5_INCLUDE_DIRS)
     set_source_files_properties(${COMMON_TEST_DIR}/hdf5_eigen.cpp PROPERTIES COMPILE_DEFINITIONS HDF5_DEFINITIONS)
 endif()
 
 add_catch_test(NAME common SOURCES
     ${SPECMICP_COMMON_TEST_FILES}
     LINK_LIBRARIES specmicp_common_static
 )
 
 # MiCPSolver
 # ----------
 
 set(MICPSOLVER_TEST_DIR micpsolver)
 
 add_catch_test(NAME micpsolver SOURCES
     ${MICPSOLVER_TEST_DIR}/test_micpsolver.cpp
     ${MICPSOLVER_TEST_DIR}/condition_number.cpp
     ${MICPSOLVER_TEST_DIR}/ncp_functions.cpp
     ${MICPSOLVER_TEST_DIR}/micpsolver.cpp
 
     LINK_LIBRARIES specmicp_common_static
 )
 
 # ODEInt
 # ----------
 
 set(ODEINT_TEST_DIR odeint)
 
 add_catch_test(NAME odeint SOURCES
     ${ODEINT_TEST_DIR}/test_odeint.cpp
     ${ODEINT_TEST_DIR}/butchertableau.cpp
     ${ODEINT_TEST_DIR}/embeddedrungekutta.cpp
 )
 
 # Database
 # --------
 
 add_subdirectory(database)
 
 
 # Specmicp
 # --------
 
 add_subdirectory( specmicp )
 
 # Reactmicp
 # ---------
 
 add_subdirectory( reactmicp )
diff --git a/tests/common/sparse_solvers.cpp b/tests/common/sparse_solvers.cpp
new file mode 100644
index 0000000..8d4c068
--- /dev/null
+++ b/tests/common/sparse_solvers.cpp
@@ -0,0 +1,299 @@
+#include "catch.hpp"
+
+#include "utils/sparse_solvers/sparse_solver.hpp"
+#include <vector>
+
+using namespace specmicp;
+
+constexpr int size = 10;
+
+using MatrixT = Eigen::SparseMatrix<scalar_t>;
+
+// simple diffusion problem
+
+// BC :
+// left : fixed value
+// right : no flow
+
+void fill_residuals(Vector& x, Vector& residuals)
+{
+
+    residuals.resize(size);
+
+    for (int i=0; i<size-1; ++i)
+    {
+        residuals(i) = (-2 *x(i+1) + x(i) + x(i+2));
+    }
+    residuals(9) = x(9)-x(10);
+}
+
+void update_solution(Vector& xp, const Vector& x,const Vector solution)
+{
+    for (int i=0; i<size; ++i)
+    {
+        xp(i+1) = x(i+1)+solution(i);
+    }
+}
+
+void fill_jacobian(MatrixT& jacobian)
+{
+    using Triplet = Eigen::Triplet<scalar_t, index_t>;
+    std::vector<Triplet> list_triplets;
+
+    list_triplets.emplace_back(Triplet(0, 0, -2));
+    list_triplets.emplace_back(Triplet(9, 9, -1));
+    for (auto i=1; i<size-1; ++i)
+    {
+        list_triplets.emplace_back(Triplet(i, i, -2));
+    }
+    for (auto i=0; i<size-1; ++i)
+    {
+        list_triplets.emplace_back(Triplet(i, i+1, 1));
+        list_triplets.emplace_back(Triplet(i+1, i, 1));
+    }
+
+    jacobian = MatrixT(size, size);
+    jacobian.setFromTriplets(list_triplets.begin(), list_triplets.end());
+
+}
+
+TEST_CASE("Sparse solvers", "[Eigen],[Sparse],[Solver]") {
+
+
+    Vector x(size+1);
+    x << 1, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2;
+
+    Vector scaling(size);
+    scaling.setConstant(1.0e3);
+
+    SECTION("LU solver", "[LU]") {
+        Vector residuals;
+        fill_residuals(x, residuals);
+
+        MatrixT jacobian;
+        fill_jacobian(jacobian);
+
+        auto lu_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
+                    sparse_solvers::SparseSolver::SparseLU);
+
+        lu_solver->analyse_pattern(jacobian);
+
+        auto retcode = lu_solver->decompose(jacobian);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector solution(size);
+        retcode = lu_solver->solve(residuals, solution);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector xp(size+1);
+        xp(0) = 1;
+        update_solution(xp, x, solution);
+        fill_residuals(xp, residuals);
+
+        REQUIRE(residuals.norm() < 1e-10);
+    }
+
+    SECTION("LU solver with scaling", "[LU],[scaling]") {
+        Vector residuals;
+        fill_residuals(x, residuals);
+
+        MatrixT jacobian;
+        fill_jacobian(jacobian);
+        jacobian = jacobian*scaling.asDiagonal();
+        jacobian.makeCompressed();
+
+        auto lu_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
+                    sparse_solvers::SparseSolver::SparseLU);
+
+        lu_solver->analyse_pattern(jacobian);
+
+        auto retcode = lu_solver->decompose(jacobian);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector solution(size);
+        retcode = lu_solver->solve_scaling(residuals, scaling, solution);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector xp(size+1);
+        xp(0) = 1;
+        update_solution(xp, x, solution);
+        fill_residuals(xp, residuals);
+
+        REQUIRE(residuals.norm() < 1e-10);
+    }
+
+
+    SECTION("QR solver", "[QR]") {
+        Vector residuals;
+        fill_residuals(x, residuals);
+
+        MatrixT jacobian;
+        fill_jacobian(jacobian);
+
+        auto qr_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
+                    sparse_solvers::SparseSolver::SparseQR);
+
+        qr_solver->analyse_pattern(jacobian);
+
+        auto retcode = qr_solver->decompose(jacobian);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector solution(size);
+        retcode = qr_solver->solve(residuals, solution);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector xp(size+1);
+        xp(0) = 1;
+        update_solution(xp, x, solution);
+        fill_residuals(xp, residuals);
+
+        REQUIRE(residuals.norm() < 1e-10);
+    }
+
+
+    SECTION("QR solver with scaling", "[QR],[scaling]") {
+        Vector residuals;
+        fill_residuals(x, residuals);
+
+        MatrixT jacobian;
+        fill_jacobian(jacobian);
+        jacobian = jacobian*scaling.asDiagonal();
+        jacobian.makeCompressed();
+
+        auto qr_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
+                    sparse_solvers::SparseSolver::SparseQR);
+
+        qr_solver->analyse_pattern(jacobian);
+
+        auto retcode = qr_solver->decompose(jacobian);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector solution(size);
+        retcode = qr_solver->solve_scaling(residuals, scaling, solution);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector xp(size+1);
+        xp(0) = 1;
+        update_solution(xp, x, solution);
+        fill_residuals(xp, residuals);
+
+        REQUIRE(residuals.norm() < 1e-10);
+    }
+
+    SECTION("BiCGSTAB solver", "[BiCGSTAB]") {
+        Vector residuals;
+        fill_residuals(x, residuals);
+
+        MatrixT jacobian;
+        fill_jacobian(jacobian);
+
+        auto bicgstab_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
+                    sparse_solvers::SparseSolver::BiCGSTAB);
+
+        bicgstab_solver->analyse_pattern(jacobian);
+
+        auto retcode = bicgstab_solver->decompose(jacobian);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector solution(size);
+        retcode = bicgstab_solver->solve(residuals, solution);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector xp(size+1);
+        xp(0) = 1;
+        update_solution(xp, x, solution);
+        fill_residuals(xp, residuals);
+
+        REQUIRE(residuals.norm() < 1e-10);
+    }
+
+    SECTION("BiCGSTAB solver with scaling", "[BiCGSTAB], [scaling]") {
+        Vector residuals;
+        fill_residuals(x, residuals);
+
+        MatrixT jacobian;
+        fill_jacobian(jacobian);
+        jacobian = jacobian*scaling.asDiagonal();
+        jacobian.makeCompressed();
+
+
+        auto bicgstab_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
+                    sparse_solvers::SparseSolver::BiCGSTAB);
+
+        bicgstab_solver->analyse_pattern(jacobian);
+
+        auto retcode = bicgstab_solver->decompose(jacobian);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector solution(size);
+        retcode = bicgstab_solver->solve_scaling(residuals, scaling, solution);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector xp(size+1);
+        xp(0) = 1;
+        update_solution(xp, x, solution);
+        fill_residuals(xp, residuals);
+
+        REQUIRE(residuals.norm() < 1e-10);
+    }
+
+#ifdef EIGEN_UNSUPPORTED_FOUND
+    SECTION("GMRES solver", "[GMRES]") {
+
+        Vector residuals;
+        fill_residuals(x, residuals);
+
+        MatrixT jacobian;
+        fill_jacobian(jacobian);
+
+        auto gmres_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
+                    sparse_solvers::SparseSolver::GMRES);
+
+        gmres_solver->analyse_pattern(jacobian);
+
+        auto retcode = gmres_solver->decompose(jacobian);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector solution(size);
+        retcode = gmres_solver->solve(residuals, solution);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector xp(size+1);
+        xp(0) = 1;
+        update_solution(xp, x, solution);
+        fill_residuals(xp, residuals);
+
+        REQUIRE(residuals.norm() < 1e-10);
+    }
+
+    SECTION("GMRES solver with scaling", "[GMRES],[scaling]") {
+
+        Vector residuals;
+        fill_residuals(x, residuals);
+
+        MatrixT jacobian;
+        fill_jacobian(jacobian);
+        jacobian = jacobian*scaling.asDiagonal();
+        jacobian.makeCompressed();
+
+        auto gmres_solver = sparse_solvers::get_sparse_solver<MatrixT, Vector, Vector>(
+                    sparse_solvers::SparseSolver::GMRES);
+
+        gmres_solver->analyse_pattern(jacobian);
+
+        auto retcode = gmres_solver->decompose(jacobian);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector solution(size);
+        retcode = gmres_solver->solve_scaling(residuals, scaling, solution);
+        REQUIRE(retcode == sparse_solvers::SparseSolverReturnCode::Success);
+
+        Vector xp(size+1);
+        xp(0) = 1;
+        update_solution(xp, x, solution);
+        fill_residuals(xp, residuals);
+
+        REQUIRE(residuals.norm() < 1e-10);
+    }
+#endif
+}