diff --git a/src/utils/sparse_solvers/sparse_bicgstab.hpp b/src/utils/sparse_solvers/sparse_bicgstab.hpp
index 585cd77..993ff64 100644
--- a/src/utils/sparse_solvers/sparse_bicgstab.hpp
+++ b/src/utils/sparse_solvers/sparse_bicgstab.hpp
@@ -1,64 +1,69 @@
 #ifndef SPECMICP_SPARSESOLVERS_SPARSEBICGSTAB_HPP
 #define SPECMICP_SPARSESOLVERS_SPARSEBICGSTAB_HPP
 
 #include "sparse_solver_structs.hpp"
 #include "sparse_solver_base.hpp"
 
 #include <Eigen/IterativeLinearSolvers>
 
 namespace specmicp {
 namespace sparse_solvers {
 
 template <class MatrixT, class DerivedR, class DerivedS>
 class SparseSolverBiCGSTAB: public SparseSolverBase<MatrixT, DerivedR, DerivedS>
 {
     using SolverT = Eigen::BiCGSTAB<MatrixT, Eigen::IncompleteLUT<typename MatrixT::Scalar>>;
 
 public:
+    void analyse_pattern(MatrixT& jacobian)
+    {
+        m_solver.analyzePattern(jacobian);
+    }
+
     //! \brief Decompose the jacboian
     SparseSolverReturnCode decompose(MatrixT& jacobian)
     {
         m_solver.compute(jacobian);
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedDecomposition;
         }
         return SparseSolverReturnCode::Success;
     }
 
     //! \brief Solve the problem
     SparseSolverReturnCode solve(
             const DerivedR& residuals,
             DerivedS& solution
             )
     {
         solution = m_solver.solve(-residuals);
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedSystemSolving;
         }
         return SparseSolverReturnCode::Success;
     }
 
     SparseSolverReturnCode solve_scaling(
             const DerivedR& residuals,
             const DerivedR& scaling,
             DerivedS& solution
             )
     {
         solution = m_solver.solve(scaling.asDiagonal()*(-residuals));
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedSystemSolving;
         }
         return SparseSolverReturnCode::Success;
     }
 
 private:
      SolverT m_solver;
 };
 
 } // end namespace sparse_solvers
 } // end namespace specmicp
 
 #endif //SPECMICP_SPARSESOLVERS_SPARSEBICGSTAB_HPP
diff --git a/src/utils/sparse_solvers/sparse_gmres.hpp b/src/utils/sparse_solvers/sparse_gmres.hpp
index a2064e3..d2a7853 100644
--- a/src/utils/sparse_solvers/sparse_gmres.hpp
+++ b/src/utils/sparse_solvers/sparse_gmres.hpp
@@ -1,67 +1,71 @@
 #ifndef SPECMICP_SPARSESOLVERS_SPARSEGMRES_HPP
 #define SPECMICP_SPARSESOLVERS_SPARSEGMRES_HPP
 
 #include "sparse_solver_structs.hpp"
 #include "sparse_solver_base.hpp"
 
 // needed for eigen 3.2.2
 #include <iostream>
 
 #include <Eigen/IterativeSolvers>
 
 namespace specmicp {
 namespace sparse_solvers {
 
 template <class MatrixT, class DerivedR, class DerivedS>
 class SparseSolverGMRES: public SparseSolverBase<MatrixT, DerivedR, DerivedS>
 {
     using SolverT = Eigen::GMRES<MatrixT, Eigen::IncompleteLUT<typename MatrixT::Scalar>> ;
 
 public:
+    void analyse_pattern(MatrixT& jacobian)
+    {
+    }
+
     //! \brief Decompose the jacboian
     SparseSolverReturnCode decompose(MatrixT& jacobian)
     {
         m_solver.compute(jacobian);
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedDecomposition;
         }
         return SparseSolverReturnCode::Success;
     }
 
     //! \brief Solve the problem
     SparseSolverReturnCode solve(
             const DerivedR& residuals,
             DerivedS& solution
             )
     {
         solution = m_solver.solve(-residuals);
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedSystemSolving;
         }
         return SparseSolverReturnCode::Success;
     }
 
     SparseSolverReturnCode solve_scaling(
             const DerivedR& residuals,
             const DerivedR& scaling,
             DerivedS& solution
             )
     {
         solution = m_solver.solve(scaling.asDiagonal()*(-residuals));
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedSystemSolving;
         }
         return SparseSolverReturnCode::Success;
     }
 
 private:
      SolverT m_solver;
 };
 
 } // end namespace sparse_solvers
 } // end namespace specmicp
 
 #endif //SPECMICP_SPARSESOLVERS_SPARSEGMRES_HPP
diff --git a/src/utils/sparse_solvers/sparse_lu.hpp b/src/utils/sparse_solvers/sparse_lu.hpp
index e59e358..0d5fb76 100644
--- a/src/utils/sparse_solvers/sparse_lu.hpp
+++ b/src/utils/sparse_solvers/sparse_lu.hpp
@@ -1,64 +1,69 @@
 #ifndef SPECMICP_SPARSESOLVERS_SPARSELU_HPP
 #define SPECMICP_SPARSESOLVERS_SPARSELU_HPP
 
 #include "sparse_solver_structs.hpp"
 #include "sparse_solver_base.hpp"
 
 #include <Eigen/SparseLU>
 
 namespace specmicp {
 namespace sparse_solvers {
 
 template <class MatrixT, class DerivedR, class DerivedS>
 class SparseSolverLU: public SparseSolverBase<MatrixT, DerivedR, DerivedS>
 {
-    using SolverT = Eigen::SparseLU<MatrixT, Eigen::COLAMDOrdering<int>>;
+    using SolverT = Eigen::SparseLU<MatrixT, Eigen::COLAMDOrdering<typename MatrixT::Index>>;
 
 public:
+    void analyse_pattern(MatrixT& jacobian)
+    {
+        m_solver.analyzePattern(jacobian);
+    }
+
     //! \brief Decompose the jacboian
     SparseSolverReturnCode decompose(MatrixT& jacobian)
     {
-        m_solver.compute(jacobian);
+        m_solver.factorize(jacobian);
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedDecomposition;
         }
         return SparseSolverReturnCode::Success;
     }
 
     //! \brief Solve the problem
     SparseSolverReturnCode solve(
             const DerivedR& residuals,
             DerivedS& solution
             )
     {
         solution = m_solver.solve(-residuals);
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedSystemSolving;
         }
         return SparseSolverReturnCode::Success;
     }
 
     SparseSolverReturnCode solve_scaling(
             const DerivedR& residuals,
             const DerivedR& scaling,
             DerivedS& solution
             )
     {
         solution = m_solver.solve(scaling.asDiagonal()*(-residuals));
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedSystemSolving;
         }
         return SparseSolverReturnCode::Success;
     }
 
 private:
      SolverT m_solver;
 };
 
 } // end namespace sparse_solvers
 } // end namespace specmicp
 
 #endif //SPECMICP_SPARSESOLVERS_SPARSELU_HPP
diff --git a/src/utils/sparse_solvers/sparse_qr.hpp b/src/utils/sparse_solvers/sparse_qr.hpp
index 470f67d..0945bf4 100644
--- a/src/utils/sparse_solvers/sparse_qr.hpp
+++ b/src/utils/sparse_solvers/sparse_qr.hpp
@@ -1,64 +1,69 @@
 #ifndef SPECMICP_SPARSESOLVERS_SPARSEQR_HPP
 #define SPECMICP_SPARSESOLVERS_SPARSEQR_HPP
 
 #include "sparse_solver_structs.hpp"
 #include "sparse_solver_base.hpp"
 
 #include <Eigen/SparseQR>
 
 namespace specmicp {
 namespace sparse_solvers {
 
 template <class MatrixT, class DerivedR, class DerivedS>
 class SparseSolverQR: public SparseSolverBase<MatrixT, DerivedR, DerivedS>
 {
     using SolverT = Eigen::SparseQR<MatrixT, Eigen::COLAMDOrdering<int>>;
 
 public:
+    void analyse_pattern(MatrixT& jacobian)
+    {
+        m_solver.analyzePattern(jacobian);
+    }
+
     //! \brief Decompose the jacboian
     SparseSolverReturnCode decompose(MatrixT& jacobian)
     {
-        m_solver = SolverT(jacobian);
+        m_solver.factorize(jacobian);
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedDecomposition;
         }
         return SparseSolverReturnCode::Success;
     }
 
     //! \brief Solve the problem
     SparseSolverReturnCode solve(
             const DerivedR& residuals,
             DerivedS& solution
             )
     {
         solution = m_solver.solve(-residuals);
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedSystemSolving;
         }
         return SparseSolverReturnCode::Success;
     }
 
     SparseSolverReturnCode solve_scaling(
             const DerivedR& residuals,
             const DerivedR& scaling,
             DerivedS& solution
             )
     {
         solution = m_solver.solve(scaling.asDiagonal()*(-residuals));
         if (m_solver.info() != Eigen::Success)
         {
             return SparseSolverReturnCode::FailedSystemSolving;
         }
         return SparseSolverReturnCode::Success;
     }
 
 private:
      SolverT m_solver;
 };
 
 } // end namespace sparse_solvers
 } // end namespace specmicp
 
 #endif //SPECMICP_SPARSESOLVERS_SPARSEQR_HPP
diff --git a/src/utils/sparse_solvers/sparse_solver_base.hpp b/src/utils/sparse_solvers/sparse_solver_base.hpp
index 58e052a..338e451 100644
--- a/src/utils/sparse_solvers/sparse_solver_base.hpp
+++ b/src/utils/sparse_solvers/sparse_solver_base.hpp
@@ -1,27 +1,28 @@
 #ifndef SPECMICP_SPARSESOLVERS_SPARSESOLVERSBASE_HPP
 #define SPECMICP_SPARSESOLVERS_SPARSESOLVERSBASE_HPP
 
 namespace specmicp {
 namespace sparse_solvers {
 
 enum class SparseSolverReturnCode;
 
 //! \brief Abstract Base Class for a sparse solver
 template <class MatrixT, class DerivedR, class DerivedS>
 class SparseSolverBase
 {
 public:
     virtual ~SparseSolverBase() {}
+    virtual void analyse_pattern(MatrixT& jacobian) = 0;
     //! \brief Decompose the jacboian
     virtual SparseSolverReturnCode decompose(MatrixT& jacobian) = 0;
     //! \brief Solve the problem
     virtual SparseSolverReturnCode solve(const DerivedR& residuals, DerivedS& solution) = 0;
     //! \brief Solve the problem
     virtual SparseSolverReturnCode solve_scaling(const DerivedR& residuals, const DerivedR& scaling, DerivedS& solution) = 0;
 };
 
 } // end namespace sparse_solvers
 } // end namespace specmicp
 
 
 #endif //SPECMICP_SPARSESOLVERS_SPARSESOLVERSBASE_HPP