diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2e62d74..16c3813 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -70,6 +70,7 @@ endif()
 option(POLYSOLVE_WITH_SANITIZERS    "Enable sanitizers in compilation targets"          OFF)
 
 # Polysolve options for enabling/disabling optional libraries
+option(POLYSOLVE_WITH_SYMPILER      "Enable sympiler library"                       ${UNIX})
 option(POLYSOLVE_WITH_ACCELERATE    "Enable Apple Accelerate" ${POLYSOLVE_ON_APPLE_SILICON})
 option(POLYSOLVE_WITH_CHOLMOD       "Enable Cholmod library"                             ON)
 option(POLYSOLVE_WITH_UMFPACK       "Enable UmfPack library"                             ON)
@@ -195,7 +196,14 @@ endif()
 include(json)
 target_link_libraries(polysolve PUBLIC nlohmann_json::nlohmann_json)
 
-# Accelerate solvers
+# Polysolve solver
+if (POLYSOLVE_WITH_SYMPILER)
+    include(sympiler)
+    target_link_libraries(polysolve PUBLIC sympiler::sym_sparse_blas)
+    target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_SYMPILER)
+endif()
+
+# Accelerate solver
 if(POLYSOLVE_WITH_ACCELERATE)
     target_compile_definitions(polysolve PUBLIC -DPOLYSOLVE_WITH_ACCELERATE)
 endif()
diff --git a/cmake/recipes/CPM.cmake b/cmake/recipes/CPM.cmake
index a3086b7..c8266bd 100644
--- a/cmake/recipes/CPM.cmake
+++ b/cmake/recipes/CPM.cmake
@@ -1,3 +1,6 @@
+# Source: https://github.com/cpm-cmake/CPM.cmake/blob/master/cmake/get_cpm.cmake
+# SPDX-License-Identifier: MIT
+
 set(CPM_DOWNLOAD_VERSION 0.38.1)
 
 if(CPM_SOURCE_CACHE)
diff --git a/cmake/recipes/sympiler.cmake b/cmake/recipes/sympiler.cmake
new file mode 100644
index 0000000..fedad67
--- /dev/null
+++ b/cmake/recipes/sympiler.cmake
@@ -0,0 +1,47 @@
+#
+# Copyright 2020 Adobe. All rights reserved.
+# This file is licensed to you under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License. You may obtain a copy
+# of the License at http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software distributed under
+# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
+# OF ANY KIND, either express or implied. See the License for the specific language
+# governing permissions and limitations under the License.
+#
+if(TARGET sympiler::sym_sparse_blas)
+    return()
+endif()
+
+message(STATUS "Third-party (external): creating target 'sympiler::sym_sparse_blas'")
+
+if("${CMAKE_SYSTEM_PROCESSOR}" MATCHES "arm64" OR "${CMAKE_OSX_ARCHITECTURES}" MATCHES "arm64")
+    # apple M1
+    set(SYMPILER_BLAS_BACKEND "Apple" CACHE STRING "BLAS implementation for SYMPILER to use")
+else()
+    # windows, linux, apple intel
+    set(SYMPILER_BLAS_BACKEND "MKL" CACHE STRING "BLAS implementation for SYMPILER to use")
+endif()
+
+if(SYMPILER_BLAS_BACKEND STREQUAL "MKL")
+    message("INCLUDING MKL")
+    include(mkl)
+    if(NOT TARGET MKL::MKL)
+        add_library(MKL::MKL ALIAS mkl::mkl)
+    endif()
+endif()
+
+include(CPM)
+CPMAddPackage(
+    NAME sympiler
+    GITHUB_REPOSITORY ryansynk/sympiler
+    GIT_TAG 2fe1ef1d72b551acedd848d64b51fa5e90251804
+)
+
+add_library(sympiler::sym_sparse_blas ALIAS sym_sparse_blas)
+
+include(GNUInstallDirs)
+target_include_directories(sym_sparse_blas INTERFACE
+    "$<BUILD_INTERFACE:${sympiler_SOURCE_DIR}>/include"
+    "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
+)
diff --git a/src/polysolve/CMakeLists.txt b/src/polysolve/CMakeLists.txt
index 831fab2..c514662 100644
--- a/src/polysolve/CMakeLists.txt
+++ b/src/polysolve/CMakeLists.txt
@@ -14,6 +14,8 @@ set(SOURCES
     LinearSolverHypre.hpp
     LinearSolverPardiso.cpp
     LinearSolverPardiso.hpp
+    LinearSolverSympiler.cpp
+    LinearSolverSympiler.hpp
     SaddlePointSolver.cpp
     SaddlePointSolver.hpp
 )
diff --git a/src/polysolve/LinearSolver.cpp b/src/polysolve/LinearSolver.cpp
index d9f0a47..3cd55f5 100644
--- a/src/polysolve/LinearSolver.cpp
+++ b/src/polysolve/LinearSolver.cpp
@@ -5,6 +5,9 @@
 
 // -----------------------------------------------------------------------------
 #include <Eigen/Sparse>
+#ifdef POLYSOLVE_WITH_SYMPILER
+#include <polysolve/LinearSolverSympiler.hpp>
+#endif
 #ifdef POLYSOLVE_WITH_ACCELERATE
 #include <Eigen/AccelerateSupport>
 #endif
@@ -191,6 +194,12 @@ namespace polysolve
         else if (solver == "Eigen::SparseLU")
         {
             RETURN_DIRECT_SOLVER_PTR(SparseLU, "Eigen::SparseLU");
+#ifdef POLYSOLVE_WITH_SYMPILER
+        }
+        else if (solver == "Sympiler")
+        {
+            return std::make_unique<LinearSolverSympiler>();
+#endif
 #ifdef POLYSOLVE_WITH_ACCELERATE
         }
         else if (solver == "Eigen::AccelerateLLT")
@@ -363,6 +372,9 @@ namespace polysolve
         return {{
             "Eigen::SimplicialLDLT",
             "Eigen::SparseLU",
+#ifdef POLYSOLVE_WITH_SYMPILER
+            "Sympiler",
+#endif
 #ifdef POLYSOLVE_WITH_ACCELERATE
             "Eigen::AccelerateLLT",
             "Eigen::AccelerateLDLT",
diff --git a/src/polysolve/LinearSolverSympiler.cpp b/src/polysolve/LinearSolverSympiler.cpp
new file mode 100644
index 0000000..54fa327
--- /dev/null
+++ b/src/polysolve/LinearSolverSympiler.cpp
@@ -0,0 +1,95 @@
+#ifdef POLYSOLVE_WITH_SYMPILER
+
+#include <polysolve/LinearSolverSympiler.hpp>
+#include <sympiler/parsy/cholesky_solver.h>
+#include <Eigen/Core>
+
+namespace polysolve
+{
+class LinearSolverSympilerImpl
+{
+public:
+    LinearSolverSympilerImpl()
+        : m_solver(nullptr), m_A_csc(std::make_unique<sym_lib::parsy::CSC>()) {}
+    ~LinearSolverSympilerImpl() = default;
+
+    std::unique_ptr<sym_lib::parsy::SolverSettings> m_solver;
+    std::unique_ptr<sym_lib::parsy::CSC> m_A_csc;
+    polysolve::StiffnessMatrix m_A_copy;
+};
+
+    LinearSolverSympiler::LinearSolverSympiler() : m_pImpl(new LinearSolverSympilerImpl) {}
+    LinearSolverSympiler::~LinearSolverSympiler() = default;
+
+    void LinearSolverSympiler::setSolverMode(int _solver_mode)
+    {
+        solver_mode = _solver_mode;
+    }
+
+    void LinearSolverSympiler::updateCSC()
+    {
+        m_pImpl->m_A_csc->nzmax = m_pImpl->m_A_copy.nonZeros();
+        m_pImpl->m_A_csc->ncol = m_pImpl->m_A_csc->nrow = m_pImpl->m_A_copy.rows();
+        m_pImpl->m_A_csc->p = m_pImpl->m_A_copy.outerIndexPtr();
+        m_pImpl->m_A_csc->i = m_pImpl->m_A_copy.innerIndexPtr();
+
+        m_pImpl->m_A_csc->x = m_pImpl->m_A_copy.valuePtr();
+        m_pImpl->m_A_csc->stype = -1;
+        m_pImpl->m_A_csc->packed = 1;
+        m_pImpl->m_A_csc->sorted = 1;
+        m_pImpl->m_A_csc->nz = nullptr;
+    }
+
+    void LinearSolverSympiler::setParameters(const json &params)
+    {
+        if (params.contains("Sympiler"))
+        {
+            if (params["Sympiler"].contains("solver_mode"))
+            {
+                setSolverMode(params["Sympiler"]["mtype"].get<int>());
+            }
+        }
+    }
+
+    void LinearSolverSympiler::getInfo(json &params) const
+    {
+        if (factorize_status == 1)
+        {
+            params["factorize_info"] = "Success";
+        }
+        else
+        {
+            params["factorize_info"] = "Failure";
+        }
+    }
+
+    void LinearSolverSympiler::analyzePattern(const StiffnessMatrix &A, const int precond_num)
+    {
+        m_pImpl->m_A_copy = StiffnessMatrix(A);
+        updateCSC();
+        m_pImpl->m_solver = std::make_unique<sym_lib::parsy::SolverSettings>(m_pImpl->m_A_csc.get());
+        m_pImpl->m_solver->symbolic_analysis();
+    }
+
+    // Factorize system matrix
+    void LinearSolverSympiler::factorize(const StiffnessMatrix &A)
+    {
+        // Only copy when matrix changes
+        if (!m_pImpl->m_A_copy.isApprox(A))
+        {
+            m_pImpl->m_A_copy = StiffnessMatrix(A);
+            updateCSC();
+        }
+        factorize_status = m_pImpl->m_solver->numerical_factorization(m_pImpl->m_A_csc.get());
+    }
+
+    // Solve the linear system Ax = b
+    void LinearSolverSympiler::solve(const Ref<const Eigen::VectorXd> b, Ref<Eigen::VectorXd> x)
+    {
+        double *x_ptr = m_pImpl->m_solver->solve_only(b.data());
+        x = Eigen::Map<Eigen::VectorXd>(x_ptr, x.rows(), x.cols());
+    }
+
+} // namespace polysolve
+
+#endif
\ No newline at end of file
diff --git a/src/polysolve/LinearSolverSympiler.hpp b/src/polysolve/LinearSolverSympiler.hpp
new file mode 100644
index 0000000..6f1b230
--- /dev/null
+++ b/src/polysolve/LinearSolverSympiler.hpp
@@ -0,0 +1,56 @@
+#ifdef POLYSOLVE_WITH_SYMPILER
+
+#include <polysolve/LinearSolver.hpp>
+
+namespace polysolve
+{
+    class LinearSolverSympilerImpl;
+    class LinearSolverSympiler : public LinearSolver
+    {
+    public:
+        LinearSolverSympiler();
+        ~LinearSolverSympiler();
+
+    private:
+        POLYSOLVE_DELETE_MOVE_COPY(LinearSolverSympiler)
+
+    protected:
+        void setSolverMode(int _solver_mode);
+        void updateCSC();
+
+    public:
+        //////////////////////
+        // Public interface //
+        //////////////////////
+
+        // Set solver parameters
+        virtual void setParameters(const json &params) override;
+
+        // Retrieve status (success/failure) of factorization
+        virtual void getInfo(json &params) const override;
+
+        // Analyze sparsity pattern
+        virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override;
+
+        // Factorize system matrix
+        virtual void factorize(const StiffnessMatrix &A) override;
+
+        // Solve the linear system Ax = b
+        virtual void solve(const Ref<const VectorXd> b, Ref<VectorXd> x) override;
+
+        // Name of the solver type (for debugging purposes)
+        virtual std::string name() const override { return "Sympiler"; }
+
+    private:
+        std::unique_ptr<LinearSolverSympilerImpl> m_pImpl;
+
+    protected:
+        ////////////////////
+        // Sympiler stuff //
+        ////////////////////
+        int solver_mode = 0; // 0 is normal solve, 1 is row/col addition
+        int factorize_status = -1;
+    };
+} // namespace polysolve
+
+#endif
\ No newline at end of file