diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
index b4c59bb0f..ade2717b0 100644
--- a/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.cpp
@@ -15,12 +15,12 @@
 namespace storm {
     namespace modelchecker {
         template<storm::dd::DdType DdType, class ValueType>
-        HybridCtmcCslModelChecker<DdType, ValueType>::HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType, ValueType> const& model) : SymbolicPropositionalModelChecker<DdType, ValueType>(model), linearEquationSolverFactory(std::make_unique<storm::utility::solver::GeneralLinearEquationSolverFactory<ValueType>>()) {
+        HybridCtmcCslModelChecker<DdType, ValueType>::HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType, ValueType> const& model) : SymbolicPropositionalModelChecker<DdType, ValueType>(model), linearEquationSolverFactory(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<ValueType>>()) {
             // Intentionally left empty.
         }
         
         template<storm::dd::DdType DdType, class ValueType>
-        HybridCtmcCslModelChecker<DdType, ValueType>::HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType, ValueType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+        HybridCtmcCslModelChecker<DdType, ValueType>::HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType, ValueType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
diff --git a/src/modelchecker/csl/HybridCtmcCslModelChecker.h b/src/modelchecker/csl/HybridCtmcCslModelChecker.h
index 5bfe49405..0d69f5220 100644
--- a/src/modelchecker/csl/HybridCtmcCslModelChecker.h
+++ b/src/modelchecker/csl/HybridCtmcCslModelChecker.h
@@ -5,7 +5,7 @@
 
 #include "src/models/symbolic/Ctmc.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
@@ -14,7 +14,7 @@ namespace storm {
         class HybridCtmcCslModelChecker : public SymbolicPropositionalModelChecker<DdType, ValueType> {
         public:
             explicit HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType, ValueType> const& model);
-            explicit HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
+            explicit HybridCtmcCslModelChecker(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
             
             // The implemented methods of the AbstractModelChecker interface.
             virtual bool canHandle(CheckTask<storm::logic::Formula> const& checkTask) const override;
@@ -31,7 +31,7 @@ namespace storm {
             
         private:
             // An object that is used for solving linear equations and performing matrix-vector multiplication.
-            std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
+            std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
         
     } // namespace modelchecker
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
index abf21fe01..f376e19e7 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
@@ -24,12 +24,12 @@
 namespace storm {
     namespace modelchecker {
         template <typename SparseCtmcModelType>
-        SparseCtmcCslModelChecker<SparseCtmcModelType>::SparseCtmcCslModelChecker(SparseCtmcModelType const& model) : SparsePropositionalModelChecker<SparseCtmcModelType>(model), linearEquationSolverFactory(std::make_unique<storm::utility::solver::GeneralLinearEquationSolverFactory<ValueType>>()) {
+        SparseCtmcCslModelChecker<SparseCtmcModelType>::SparseCtmcCslModelChecker(SparseCtmcModelType const& model) : SparsePropositionalModelChecker<SparseCtmcModelType>(model), linearEquationSolverFactory(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<ValueType>>()) {
             // Intentionally left empty.
         }
         
         template <typename SparseCtmcModelType>
-        SparseCtmcCslModelChecker<SparseCtmcModelType>::SparseCtmcCslModelChecker(SparseCtmcModelType const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SparsePropositionalModelChecker<SparseCtmcModelType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+        SparseCtmcCslModelChecker<SparseCtmcModelType>::SparseCtmcCslModelChecker(SparseCtmcModelType const& model, std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SparsePropositionalModelChecker<SparseCtmcModelType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
index 802841dee..070f81ffd 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
@@ -5,7 +5,7 @@
 
 #include "src/models/sparse/Ctmc.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
@@ -17,7 +17,7 @@ namespace storm {
             typedef typename SparseCtmcModelType::RewardModelType RewardModelType;
             
             explicit SparseCtmcCslModelChecker(SparseCtmcModelType const& model);
-            explicit SparseCtmcCslModelChecker(SparseCtmcModelType const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
+            explicit SparseCtmcCslModelChecker(SparseCtmcModelType const& model, std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
             
             // The implemented methods of the AbstractModelChecker interface.
             virtual bool canHandle(CheckTask<storm::logic::Formula> const& checkTask) const override;
@@ -32,7 +32,7 @@ namespace storm {
 
         private:
             // An object that is used for solving linear equations and performing matrix-vector multiplication.
-            std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
+            std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
         
     } // namespace modelchecker
diff --git a/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp b/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp
index e59523866..91581e987 100644
--- a/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp
+++ b/src/modelchecker/csl/helper/HybridCtmcCslHelper.cpp
@@ -26,7 +26,7 @@ namespace storm {
         namespace helper {
             
             template<storm::dd::DdType DdType, class ValueType>
-            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeReachabilityRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeReachabilityRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 
                 return HybridDtmcPrctlHelper<DdType, ValueType>::computeReachabilityRewards(model, computeProbabilityMatrix(model, rateMatrix, exitRateVector), rewardModel.divideStateRewardVector(exitRateVector), targetStates, qualitative, linearEquationSolverFactory);
             }
@@ -37,12 +37,12 @@ namespace storm {
             }
             
             template<storm::dd::DdType DdType, class ValueType>
-            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeUntilProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeUntilProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 return HybridDtmcPrctlHelper<DdType, ValueType>::computeUntilProbabilities(model, computeProbabilityMatrix(model, rateMatrix, exitRateVector), phiStates, psiStates, qualitative, linearEquationSolverFactory);
             }
             
             template<storm::dd::DdType DdType, class ValueType>
-            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeBoundedUntilProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, double lowerBound, double upperBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeBoundedUntilProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, double lowerBound, double upperBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 
                 // If the time bounds are [0, inf], we rather call untimed reachability.
                 if (storm::utility::isZero(lowerBound) && upperBound == storm::utility::infinity<ValueType>()) {
@@ -218,7 +218,7 @@ namespace storm {
             }
             
             template<storm::dd::DdType DdType, class ValueType>
-            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeInstantaneousRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeInstantaneousRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 
                 // Only compute the result if the model has a state-based reward model.
                 STORM_LOG_THROW(rewardModel.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
@@ -244,7 +244,7 @@ namespace storm {
             }
             
             template<storm::dd::DdType DdType, class ValueType>
-            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeCumulativeRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeCumulativeRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Only compute the result if the model has a state-based reward model.
                 STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
                 
@@ -276,7 +276,7 @@ namespace storm {
             }
             
             template<storm::dd::DdType DdType, class ValueType>
-            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeLongRunAverageProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridCtmcCslHelper<DdType, ValueType>::computeLongRunAverageProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 storm::dd::Add<DdType, ValueType> probabilityMatrix = computeProbabilityMatrix(model, rateMatrix, exitRateVector);
                 
                 // Create ODD for the translation.
diff --git a/src/modelchecker/csl/helper/HybridCtmcCslHelper.h b/src/modelchecker/csl/helper/HybridCtmcCslHelper.h
index ca7d7aa35..e2e0667c0 100644
--- a/src/modelchecker/csl/helper/HybridCtmcCslHelper.h
+++ b/src/modelchecker/csl/helper/HybridCtmcCslHelper.h
@@ -7,7 +7,7 @@
 
 #include "src/modelchecker/results/CheckResult.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
@@ -18,17 +18,17 @@ namespace storm {
             public:
                 typedef typename storm::models::symbolic::Model<DdType, ValueType>::RewardModelType RewardModelType;
 
-                static std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, double lowerBound, double upperBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, double lowerBound, double upperBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeCumulativeRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeCumulativeRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeUntilProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeUntilProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeReachabilityRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeReachabilityRewards(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeLongRunAverageProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeLongRunAverageProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
                 static std::unique_ptr<CheckResult> computeNextProbabilities(storm::models::symbolic::Ctmc<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& rateMatrix, storm::dd::Add<DdType, ValueType> const& exitRateVector, storm::dd::Bdd<DdType> const& nextStates);
 
diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
index 12bb227b9..0dddea5cb 100644
--- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
+++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
@@ -26,7 +26,7 @@ namespace storm {
     namespace modelchecker {
         namespace helper {
             template <typename ValueType>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 
                 uint_fast64_t numberOfStates = rateMatrix.getRowCount();
                 
@@ -191,7 +191,7 @@ namespace storm {
             }
 
             template <typename ValueType>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 return SparseDtmcPrctlHelper<ValueType>::computeUntilProbabilities(computeProbabilityMatrix(rateMatrix, exitRateVector), backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolverFactory);
             }
             
@@ -202,7 +202,7 @@ namespace storm {
             }
 
             template <typename ValueType>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 return SparseDtmcPrctlHelper<ValueType>::computeNextProbabilities(computeProbabilityMatrix(rateMatrix, exitRateVector), nextStates, linearEquationSolverFactory);
             }
             
@@ -235,7 +235,7 @@ namespace storm {
             
             template <typename ValueType>
             template<bool computeCumulativeReward>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 
                 ValueType lambda = timeBound * uniformizationRate;
                 
@@ -314,7 +314,7 @@ namespace storm {
             
             template <typename ValueType>
             template <typename RewardModelType>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Only compute the result if the model has a state-based reward this->getModel().
                 STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
                 
@@ -341,7 +341,7 @@ namespace storm {
             
             template <typename ValueType>
             template <typename RewardModelType>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Only compute the result if the model has a state-based reward this->getModel().
                 STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
                                 
@@ -373,7 +373,7 @@ namespace storm {
 
             template <typename ValueType>
             template <typename RewardModelType>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
 
                 storm::storage::SparseMatrix<ValueType> probabilityMatrix = computeProbabilityMatrix(rateMatrix, exitRateVector);
@@ -432,7 +432,7 @@ namespace storm {
             }
 
             template <typename ValueType>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // If there are no goal states, we avoid the computation and directly return zero.
                 uint_fast64_t numberOfStates = probabilityMatrix.getRowCount();
                 if (psiStates.empty()) {
@@ -656,7 +656,7 @@ namespace storm {
             }
 
             template <typename ValueType>
-            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeReachabilityTimes(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseCtmcCslHelper<ValueType>::computeReachabilityTimes(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Compute expected time on CTMC by reduction to DTMC with rewards.
                 storm::storage::SparseMatrix<ValueType> probabilityMatrix = computeProbabilityMatrix(rateMatrix, exitRateVector);
 
@@ -698,16 +698,16 @@ namespace storm {
             }
 
             template class SparseCtmcCslHelper<double>;
-            template std::vector<double> SparseCtmcCslHelper<double>::computeInstantaneousRewards(storm::storage::SparseMatrix<double> const& rateMatrix, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<double> const& linearEquationSolverFactory);
-            template std::vector<double> SparseCtmcCslHelper<double>::computeCumulativeRewards(storm::storage::SparseMatrix<double> const& rateMatrix, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<double> const& linearEquationSolverFactory);
-            template std::vector<double> SparseCtmcCslHelper<double>::computeReachabilityRewards(storm::storage::SparseMatrix<double> const& rateMatrix, storm::storage::SparseMatrix<double> const& backwardTransitions, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<double> const& linearEquationSolverFactory);
+            template std::vector<double> SparseCtmcCslHelper<double>::computeInstantaneousRewards(storm::storage::SparseMatrix<double> const& rateMatrix, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<double> const& linearEquationSolverFactory);
+            template std::vector<double> SparseCtmcCslHelper<double>::computeCumulativeRewards(storm::storage::SparseMatrix<double> const& rateMatrix, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<double> const& linearEquationSolverFactory);
+            template std::vector<double> SparseCtmcCslHelper<double>::computeReachabilityRewards(storm::storage::SparseMatrix<double> const& rateMatrix, storm::storage::SparseMatrix<double> const& backwardTransitions, std::vector<double> const& exitRateVector, storm::models::sparse::StandardRewardModel<double> const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<double> const& linearEquationSolverFactory);
 
-            template std::vector<storm::RationalNumber> SparseCtmcCslHelper<storm::RationalNumber>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<storm::RationalNumber> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<storm::RationalNumber> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<storm::RationalNumber> const& linearEquationSolverFactory);
+            template std::vector<storm::RationalNumber> SparseCtmcCslHelper<storm::RationalNumber>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<storm::RationalNumber> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<storm::RationalNumber> const* exitRateVector, bool qualitative, storm::solver::LinearEquationSolverFactory<storm::RationalNumber> const& linearEquationSolverFactory);
 
             
             template std::vector<storm::RationalFunction> SparseCtmcCslHelper<storm::RationalFunction>::computeUntilProbabilitiesElimination(storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, storm::storage::SparseMatrix<storm::RationalFunction> const& backwardTransitions, std::vector<storm::RationalFunction> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative);
             template std::vector<storm::RationalFunction> SparseCtmcCslHelper<storm::RationalFunction>::computeReachabilityTimesElimination(storm::storage::SparseMatrix<storm::RationalFunction> const& rateMatrix, storm::storage::SparseMatrix<storm::RationalFunction> const& backwardTransitions, std::vector<storm::RationalFunction> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, bool qualitative);
-            template std::vector<storm::RationalFunction> SparseCtmcCslHelper<storm::RationalFunction>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<storm::RationalFunction> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<storm::RationalFunction> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<storm::RationalFunction> const& linearEquationSolverFactory);
+            template std::vector<storm::RationalFunction> SparseCtmcCslHelper<storm::RationalFunction>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<storm::RationalFunction> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<storm::RationalFunction> const* exitRateVector, bool qualitative, storm::solver::LinearEquationSolverFactory<storm::RationalFunction> const& linearEquationSolverFactory);
         }
     }
 }
diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.h b/src/modelchecker/csl/helper/SparseCtmcCslHelper.h
index 2b4dd6043..8485b2501 100644
--- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.h
+++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.h
@@ -3,7 +3,7 @@
 
 #include "src/storage/BitVector.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
@@ -11,25 +11,25 @@ namespace storm {
             template <typename ValueType>
             class SparseCtmcCslHelper {
             public:
-                static std::vector<ValueType> computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
 
-                static std::vector<ValueType> computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 static std::vector<ValueType> computeUntilProbabilitiesElimination(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative);
                 
                 template <typename RewardModelType>
-                static std::vector<ValueType> computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
                 template <typename RewardModelType>
-                static std::vector<ValueType> computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, double timeBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
                 template <typename RewardModelType>
-                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& probabilityMatrix, storm::storage::BitVector const& psiStates, std::vector<ValueType> const* exitRateVector, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeReachabilityTimes(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
+                static std::vector<ValueType> computeReachabilityTimes(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
                 static std::vector<ValueType> computeReachabilityTimesElimination(storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& targetStates, bool qualitative);
 
                 /*!
@@ -58,7 +58,7 @@ namespace storm {
                  * @return The vector of transient probabilities.
                  */
                 template<bool useMixedPoissonProbabilities = false>
-                static std::vector<ValueType> computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, std::vector<ValueType> const* addVector, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
                 /*!
                  * Converts the given rate-matrix into a time-abstract probability matrix.
diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
index 75f651e8a..f0cd46154 100644
--- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
@@ -26,12 +26,12 @@
 namespace storm {
     namespace modelchecker {
         template<storm::dd::DdType DdType, typename ValueType>
-        HybridDtmcPrctlModelChecker<DdType, ValueType>::HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType, ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType, ValueType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+        HybridDtmcPrctlModelChecker<DdType, ValueType>::HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType, ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SymbolicPropositionalModelChecker<DdType, ValueType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        HybridDtmcPrctlModelChecker<DdType, ValueType>::HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType, ValueType> const& model) : SymbolicPropositionalModelChecker<DdType, ValueType>(model), linearEquationSolverFactory(std::make_unique<storm::utility::solver::GeneralLinearEquationSolverFactory<ValueType>>()) {
+        HybridDtmcPrctlModelChecker<DdType, ValueType>::HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType, ValueType> const& model) : SymbolicPropositionalModelChecker<DdType, ValueType>(model), linearEquationSolverFactory(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<ValueType>>()) {
             // Intentionally left empty.
         }
         
diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
index 3489e74e8..620aa2952 100644
--- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h
@@ -5,7 +5,7 @@
 
 #include "src/models/symbolic/Dtmc.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
@@ -14,7 +14,7 @@ namespace storm {
         class HybridDtmcPrctlModelChecker : public SymbolicPropositionalModelChecker<DdType, ValueType> {
         public:
             explicit HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType, ValueType> const& model);
-            explicit HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType, ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
+            explicit HybridDtmcPrctlModelChecker(storm::models::symbolic::Dtmc<DdType, ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
             
             // The implemented methods of the AbstractModelChecker interface.
             virtual bool canHandle(CheckTask<storm::logic::Formula> const& checkTask) const override;
@@ -32,7 +32,7 @@ namespace storm {
             
         private:
             // An object that is used for retrieving linear equation solvers.
-            std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
+            std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
         
     } // namespace modelchecker
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index 75fdaca64..92ec6d2be 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -24,12 +24,12 @@
 namespace storm {
     namespace modelchecker {
         template<typename SparseDtmcModelType>
-        SparseDtmcPrctlModelChecker<SparseDtmcModelType>::SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SparsePropositionalModelChecker<SparseDtmcModelType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
+        SparseDtmcPrctlModelChecker<SparseDtmcModelType>::SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model, std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SparsePropositionalModelChecker<SparseDtmcModelType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
         template<typename SparseDtmcModelType>
-        SparseDtmcPrctlModelChecker<SparseDtmcModelType>::SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model) : SparsePropositionalModelChecker<SparseDtmcModelType>(model), linearEquationSolverFactory(std::make_unique<storm::utility::solver::GeneralLinearEquationSolverFactory<ValueType>>()) {
+        SparseDtmcPrctlModelChecker<SparseDtmcModelType>::SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model) : SparsePropositionalModelChecker<SparseDtmcModelType>(model), linearEquationSolverFactory(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<ValueType>>()) {
             // Intentionally left empty.
         }
         
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
index 6b94f6642..d577676b0 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
@@ -16,7 +16,7 @@ namespace storm {
             typedef typename SparseDtmcModelType::RewardModelType RewardModelType;
                         
             explicit SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model);
-            explicit SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
+            explicit SparseDtmcPrctlModelChecker(SparseDtmcModelType const& model, std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
             
             // The implemented methods of the AbstractModelChecker interface.
             virtual bool canHandle(CheckTask<storm::logic::Formula> const& checkTask) const override;
@@ -34,7 +34,7 @@ namespace storm {
 
         private:
             // An object that is used for retrieving linear equation solvers.
-            std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
+            std::unique_ptr<storm::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
         
     } // namespace modelchecker
diff --git a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp
index 869242db4..73bd63088 100644
--- a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp
@@ -24,7 +24,7 @@ namespace storm {
         namespace helper {
 
             template<storm::dd::DdType DdType, typename ValueType>
-            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeUntilProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeUntilProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
                 // probability 0 and 1 of satisfying the until-formula.
                 std::pair<storm::dd::Bdd<DdType>, storm::dd::Bdd<DdType>> statesWithProbability01 = storm::utility::graph::performProb01(model, transitionMatrix, phiStates, psiStates);
@@ -83,7 +83,7 @@ namespace storm {
             }
 
             template<storm::dd::DdType DdType, typename ValueType>
-            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeGloballyProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeGloballyProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 std::unique_ptr<CheckResult> result = computeUntilProbabilities(model, transitionMatrix, model.getReachableStates(), !psiStates && model.getReachableStates(), qualitative, linearEquationSolverFactory);
                 result->asQuantitativeCheckResult().oneMinus();
                 return result;
@@ -96,7 +96,7 @@ namespace storm {
             }
 
             template<storm::dd::DdType DdType, typename ValueType>
-            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeBoundedUntilProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeBoundedUntilProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // We need to identify the states which have to be taken out of the matrix, i.e. all states that have
                 // probability 0 or 1 of satisfying the until-formula.
                 storm::dd::Bdd<DdType> statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(model, transitionMatrix.notZero(), phiStates, psiStates, stepBound);
@@ -141,7 +141,7 @@ namespace storm {
 
             
             template<storm::dd::DdType DdType, typename ValueType>
-            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeInstantaneousRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeInstantaneousRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Only compute the result if the model has at least one reward this->getModel().
                 STORM_LOG_THROW(rewardModel.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
                 
@@ -163,7 +163,7 @@ namespace storm {
             }
             
             template<storm::dd::DdType DdType, typename ValueType>
-            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeCumulativeRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeCumulativeRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Only compute the result if the model has at least one reward this->getModel().
                 STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
                 
@@ -189,7 +189,7 @@ namespace storm {
             }
             
             template<storm::dd::DdType DdType, typename ValueType>
-            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeReachabilityRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::unique_ptr<CheckResult> HybridDtmcPrctlHelper<DdType, ValueType>::computeReachabilityRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 
                 // Only compute the result if there is at least one reward model.
                 STORM_LOG_THROW(!rewardModel.empty(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
diff --git a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h
index 9cbfa83d3..6296802ae 100644
--- a/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h
+++ b/src/modelchecker/prctl/helper/HybridDtmcPrctlHelper.h
@@ -6,7 +6,7 @@
 #include "src/storage/dd/Add.h"
 #include "src/storage/dd/Bdd.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
@@ -20,19 +20,19 @@ namespace storm {
             public:
                 typedef typename storm::models::symbolic::Model<DdType, ValueType>::RewardModelType RewardModelType;
                 
-                static std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeBoundedUntilProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
                 static std::unique_ptr<CheckResult> computeNextProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& nextStates);
                 
-                static std::unique_ptr<CheckResult> computeUntilProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeUntilProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& phiStates, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeGloballyProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeGloballyProbabilities(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeCumulativeRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeCumulativeRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeInstantaneousRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::unique_ptr<CheckResult> computeReachabilityRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::unique_ptr<CheckResult> computeReachabilityRewards(storm::models::symbolic::Model<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, RewardModelType const& rewardModel, storm::dd::Bdd<DdType> const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             };
             
         }
diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
index 857da9a6f..7e9d03c6d 100644
--- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
@@ -17,7 +17,7 @@ namespace storm {
     namespace modelchecker {
         namespace helper {
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::zero<ValueType>());
                 
                 // If we identify the states that have probability 0 of reaching the target states, we can exclude them in the further analysis.
@@ -48,7 +48,7 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // We need to identify the states which have to be taken out of the matrix, i.e.
                 // all states that have probability 0 and 1 of satisfying the until-formula.
                 std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = storm::utility::graph::performProb01(backwardTransitions, phiStates, psiStates);
@@ -105,7 +105,7 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeGloballyProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeGloballyProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 std::vector<ValueType> result = computeUntilProbabilities(transitionMatrix, backwardTransitions, storm::storage::BitVector(transitionMatrix.getRowCount(), true), ~psiStates, qualitative, linearEquationSolverFactory);
                 for (auto& entry : result) {
                     entry = storm::utility::one<ValueType>() - entry;
@@ -114,7 +114,7 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Create the vector with which to multiply and initialize it correctly.
                 std::vector<ValueType> result(transitionMatrix.getRowCount());
                 storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
@@ -126,7 +126,7 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Initialize result to the null vector.
                 std::vector<ValueType> result(transitionMatrix.getRowCount());
                 
@@ -141,7 +141,7 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepCount, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepCount, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Only compute the result if the model has a state-based reward this->getModel().
                 STORM_LOG_THROW(rewardModel.hasStateRewards(), storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
                 
@@ -156,12 +156,12 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 return computeReachabilityRewards(transitionMatrix, backwardTransitions, [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates) { return rewardModel.getTotalRewardVector(numberOfRows, transitionMatrix, maybeStates); }, targetStates, qualitative, linearEquationSolverFactory);
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& totalStateRewardVector, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& totalStateRewardVector, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 return computeReachabilityRewards(transitionMatrix, backwardTransitions,
                                                   [&] (uint_fast64_t numberOfRows, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& maybeStates) {
                                                       std::vector<ValueType> result(numberOfRows);
@@ -172,7 +172,7 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                                 
                 // Determine which states have a reward of infinity by definition.
                 storm::storage::BitVector trueStates(transitionMatrix.getRowCount(), true);
@@ -224,12 +224,12 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 return SparseCtmcCslHelper<ValueType>::computeLongRunAverageProbabilities(transitionMatrix, psiStates, nullptr, qualitative, linearEquationSolverFactory);
             }
             
             template<typename ValueType, typename RewardModelType>
-            typename SparseDtmcPrctlHelper<ValueType, RewardModelType>::BaierTransformedModel SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeBaierTransformation(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, boost::optional<std::vector<ValueType>> const& stateRewards, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            typename SparseDtmcPrctlHelper<ValueType, RewardModelType>::BaierTransformedModel SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeBaierTransformation(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, boost::optional<std::vector<ValueType>> const& stateRewards, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 
                 BaierTransformedModel result;
                 
@@ -349,7 +349,7 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeConditionalProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeConditionalProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 
                 // Prepare result vector.
                 std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::infinity<ValueType>());
@@ -376,7 +376,7 @@ namespace storm {
             }
             
             template<typename ValueType, typename RewardModelType>
-            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeConditionalRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
+            std::vector<ValueType> SparseDtmcPrctlHelper<ValueType, RewardModelType>::computeConditionalRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) {
                 // Prepare result vector.
                 std::vector<ValueType> result(transitionMatrix.getRowCount(), storm::utility::infinity<ValueType>());
                 
diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
index 0be762298..80e4321ba 100644
--- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
+++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h
@@ -8,7 +8,7 @@
 #include "src/storage/SparseMatrix.h"
 #include "src/storage/BitVector.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
     namespace modelchecker {
@@ -19,30 +19,30 @@ namespace storm {
             template <typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
             class SparseDtmcPrctlHelper {
             public:
-                static std::vector<ValueType> computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeBoundedUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeNextProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeUntilProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeGloballyProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeGloballyProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeCumulativeRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepCount, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeInstantaneousRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepCount, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& totalStateRewardVector, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& totalStateRewardVector, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeLongRunAverageProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeConditionalProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeConditionalProbabilities(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
-                static std::vector<ValueType> computeConditionalRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeConditionalRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, RewardModelType const& rewardModel, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
             private:
-                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static std::vector<ValueType> computeReachabilityRewards(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::function<std::vector<ValueType>(uint_fast64_t, storm::storage::SparseMatrix<ValueType> const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::BitVector const& targetStates, bool qualitative, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
                 
                 struct BaierTransformedModel {
                     BaierTransformedModel() : noTargetStates(false) {
@@ -56,7 +56,7 @@ namespace storm {
                     bool noTargetStates;
                 };
                 
-                static BaierTransformedModel computeBaierTransformation(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, boost::optional<std::vector<ValueType>> const& stateRewards, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+                static BaierTransformedModel computeBaierTransformation(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::storage::BitVector const& conditionStates, boost::optional<std::vector<ValueType>> const& stateRewards, storm::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
             };
             
         }
diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp
index d7328f8f9..3d5f5a857 100644
--- a/src/solver/EigenLinearEquationSolver.cpp
+++ b/src/solver/EigenLinearEquationSolver.cpp
@@ -5,16 +5,14 @@
 #include "src/settings/SettingsManager.h"
 #include "src/settings/modules/EigenEquationSolverSettings.h"
 
+#include "src/utility/macros.h"
+#include "src/exceptions/InvalidSettingsException.h"
+
 namespace storm {
     namespace solver {
      
         template<typename ValueType>
-        EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A)), method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner) {
-            // Intentionally left empty.
-        }
-        
-        template<typename ValueType>
-        EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A)) {
+        EigenLinearEquationSolverSettings<ValueType>::EigenLinearEquationSolverSettings() {
             // Get the settings object to customize linear solving.
             storm::settings::modules::EigenEquationSolverSettings const& settings = storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>();
             
@@ -41,6 +39,65 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void EigenLinearEquationSolverSettings<ValueType>::setSolutionMethod(SolutionMethod const& method) {
+            this->method = method;
+        }
+        
+        template<typename ValueType>
+        void EigenLinearEquationSolverSettings<ValueType>::setPreconditioner(Preconditioner const& preconditioner) {
+            this->preconditioner = preconditioner;
+        }
+        
+        template<typename ValueType>
+        void EigenLinearEquationSolverSettings<ValueType>::setPrecision(ValueType precision) {
+            this->precision = precision;
+        }
+        
+        template<typename ValueType>
+        void EigenLinearEquationSolverSettings<ValueType>::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) {
+            this->maximalNumberOfIterations = maximalNumberOfIterations;
+        }
+        
+        template<typename ValueType>
+        typename EigenLinearEquationSolverSettings<ValueType>::SolutionMethod EigenLinearEquationSolverSettings<ValueType>::getSolutionMethod() const {
+            return this->method;
+        }
+        
+        template<typename ValueType>
+        typename EigenLinearEquationSolverSettings<ValueType>::Preconditioner EigenLinearEquationSolverSettings<ValueType>::getPreconditioner() const {
+            return this->preconditioner;
+        }
+        
+        template<typename ValueType>
+        ValueType EigenLinearEquationSolverSettings<ValueType>::getPrecision() const {
+            return this->precision;
+        }
+        
+        template<typename ValueType>
+        uint64_t EigenLinearEquationSolverSettings<ValueType>::getMaximalNumberOfIterations() const {
+            return this->maximalNumberOfIterations;
+        }
+        
+        EigenLinearEquationSolverSettings<storm::RationalNumber>::EigenLinearEquationSolverSettings() {
+            // Intentionally left empty.
+        }
+
+        EigenLinearEquationSolverSettings<storm::RationalFunction>::EigenLinearEquationSolverSettings() {
+            // Intentionally left empty.
+        }
+
+        template<typename ValueType>
+        EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EigenLinearEquationSolverSettings<ValueType> const& settings) : eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A)), settings(settings) {
+            // Intentionally left empty.
+        }
+
+        template<typename ValueType>
+        EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EigenLinearEquationSolverSettings<ValueType> const& settings) : settings(settings) {
+            storm::storage::SparseMatrix<ValueType> localA(std::move(A));
+            eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(localA);
+        }
+        
         template<typename ValueType>
         void EigenLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             // Translate the vectors x and b into Eigen's format.
@@ -54,25 +111,25 @@ namespace storm {
                 eigenX[index] = x[index];
             }
             
-            if (method == SolutionMethod::SparseLU) {
+            if (this->getSettings().getSolutionMethod() == EigenLinearEquationSolverSettings<ValueType>::SolutionMethod::SparseLU) {
                 Eigen::SparseLU<Eigen::SparseMatrix<ValueType>, Eigen::COLAMDOrdering<int>> solver;
-                solver.compute(*eigenA);
+                solver.compute(*this->eigenA);
                 if (solver.info() != Eigen::Success) {
                     std::cout << solver.lastErrorMessage() << std::endl;
                 }
                 solver._solve_impl(eigenB, eigenX);
             } else {
-                if (preconditioner == Preconditioner::Ilu) {
+                if (this->getSettings().getPreconditioner() == EigenLinearEquationSolverSettings<ValueType>::Preconditioner::Ilu) {
                     Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::IncompleteLUT<ValueType>> solver;
-                    solver.compute(*eigenA);
+                    solver.compute(*this->eigenA);
                     solver.solveWithGuess(eigenB, eigenX);
-                } else if (preconditioner == Preconditioner::Diagonal) {
+                } else if (this->getSettings().getPreconditioner() == EigenLinearEquationSolverSettings<ValueType>::Preconditioner::Diagonal) {
                     Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::DiagonalPreconditioner<ValueType>> solver;
-                    solver.compute(*eigenA);
+                    solver.compute(*this->eigenA);
                     solver.solveWithGuess(eigenB, eigenX);
                 } else {
                     Eigen::BiCGSTAB<Eigen::SparseMatrix<ValueType>, Eigen::IdentityPreconditioner> solver;
-                    solver.compute(*eigenA);
+                    solver.compute(*this->eigenA);
                     solver.solveWithGuess(eigenB, eigenX);
                 }
             }
@@ -100,7 +157,7 @@ namespace storm {
             
             // Perform n matrix-vector multiplications.
             for (uint64_t iteration = 0; iteration < n; ++iteration) {
-                eigenX = *eigenA * eigenX;
+                eigenX = *this->eigenA * eigenX;
                 if (eigenB != nullptr) {
                     eigenX += *eigenB;
                 }
@@ -112,12 +169,19 @@ namespace storm {
             }
         }
         
-        // Specialization form storm::RationalNumber
+        template<typename ValueType>
+        EigenLinearEquationSolverSettings<ValueType>& EigenLinearEquationSolver<ValueType>::getSettings() {
+            return settings;
+        }
         
-        EigenLinearEquationSolver<storm::RationalNumber>::EigenLinearEquationSolver(storm::storage::SparseMatrix<storm::RationalNumber> const& A, SolutionMethod method) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<storm::RationalNumber>(A)), method(method) {
-            // Intentionally left empty.
+        template<typename ValueType>
+        EigenLinearEquationSolverSettings<ValueType> const& EigenLinearEquationSolver<ValueType>::getSettings() const {
+            return settings;
         }
         
+        // Specialization form storm::RationalNumber
+        
+        template<>
         void EigenLinearEquationSolver<storm::RationalNumber>::solveEquationSystem(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>* multiplyResult) const {
             // Translate the vectors x and b into Eigen's format.
             Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1> eigenB(b.size());
@@ -130,11 +194,9 @@ namespace storm {
                 eigenX[index] = x[index];
             }
             
-            if (method == SolutionMethod::SparseLU) {
-                Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalNumber>, Eigen::COLAMDOrdering<int>> solver;
-                solver.compute(*eigenA);
-                solver._solve_impl(eigenB, eigenX);
-            }
+            Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalNumber>, Eigen::COLAMDOrdering<int>> solver;
+            solver.compute(*eigenA);
+            solver._solve_impl(eigenB, eigenX);
             
             // Translate the solution from Eigen's format into our representation.
             for (uint64_t index = 0; index < eigenX.size(); ++index) {
@@ -142,6 +204,7 @@ namespace storm {
             }
         }
         
+        template<>
         void EigenLinearEquationSolver<storm::RationalNumber>::performMatrixVectorMultiplication(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const* b, uint_fast64_t n, std::vector<storm::RationalNumber>* multiplyResult) const {
             // Translate the vectors x and b into Eigen's format.
             std::unique_ptr<Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>> eigenB;
@@ -172,10 +235,7 @@ namespace storm {
         
         // Specialization form storm::RationalFunction
         
-        EigenLinearEquationSolver<storm::RationalFunction>::EigenLinearEquationSolver(storm::storage::SparseMatrix<storm::RationalFunction> const& A, SolutionMethod method) : originalA(&A), eigenA(storm::adapters::EigenAdapter::toEigenSparseMatrix<storm::RationalFunction>(A)), method(method) {
-            // Intentionally left empty.
-        }
-        
+        template<>
         void EigenLinearEquationSolver<storm::RationalFunction>::solveEquationSystem(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b, std::vector<storm::RationalFunction>* multiplyResult) const {
             // Translate the vectors x and b into Eigen's format.
             Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1> eigenB(b.size());
@@ -188,11 +248,9 @@ namespace storm {
                 eigenX[index] = x[index];
             }
             
-            if (method == SolutionMethod::SparseLU) {
-                Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalFunction>, Eigen::COLAMDOrdering<int>> solver;
-                solver.compute(*eigenA);
-                solver._solve_impl(eigenB, eigenX);
-            }
+            Eigen::SparseLU<Eigen::SparseMatrix<storm::RationalFunction>, Eigen::COLAMDOrdering<int>> solver;
+            solver.compute(*eigenA);
+            solver._solve_impl(eigenB, eigenX);
             
             // Translate the solution from Eigen's format into our representation.
             for (uint64_t index = 0; index < eigenX.size(); ++index) {
@@ -200,6 +258,7 @@ namespace storm {
             }
         }
         
+        template<>
         void EigenLinearEquationSolver<storm::RationalFunction>::performMatrixVectorMultiplication(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const* b, uint_fast64_t n, std::vector<storm::RationalFunction>* multiplyResult) const {
             // Translate the vectors x and b into Eigen's format.
             std::unique_ptr<Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>> eigenB;
@@ -228,9 +287,37 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EigenLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+            return std::make_unique<storm::solver::EigenLinearEquationSolver<ValueType>>(matrix, settings);
+        }
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EigenLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
+            return std::make_unique<storm::solver::EigenLinearEquationSolver<ValueType>>(std::move(matrix), settings);
+        }
+        
+        template<typename ValueType>
+        EigenLinearEquationSolverSettings<ValueType>& EigenLinearEquationSolverFactory<ValueType>::getSettings() {
+            return settings;
+        }
+        
+        template<typename ValueType>
+        EigenLinearEquationSolverSettings<ValueType> const& EigenLinearEquationSolverFactory<ValueType>::getSettings() const {
+            return settings;
+        }
+        
+        template class EigenLinearEquationSolverSettings<double>;
+        template class EigenLinearEquationSolverSettings<storm::RationalNumber>;
+        template class EigenLinearEquationSolverSettings<storm::RationalFunction>;
+        
         template class EigenLinearEquationSolver<double>;
         template class EigenLinearEquationSolver<storm::RationalNumber>;
-//        template class EigenLinearEquationSolver<storm::RationalFunction>;
+        template class EigenLinearEquationSolver<storm::RationalFunction>;
+        
+        template class EigenLinearEquationSolverFactory<double>;
+        template class EigenLinearEquationSolverFactory<storm::RationalNumber>;
+        template class EigenLinearEquationSolverFactory<storm::RationalFunction>;
         
     }
 }
\ No newline at end of file
diff --git a/src/solver/EigenLinearEquationSolver.h b/src/solver/EigenLinearEquationSolver.h
index 924a2af74..6aa770d7e 100644
--- a/src/solver/EigenLinearEquationSolver.h
+++ b/src/solver/EigenLinearEquationSolver.h
@@ -7,11 +7,8 @@
 namespace storm {
     namespace solver {
         
-        /*!
-         * A class that uses the Eigen library to implement the LinearEquationSolver interface.
-         */
         template<typename ValueType>
-        class EigenLinearEquationSolver : public LinearEquationSolver<ValueType> {
+        class EigenLinearEquationSolverSettings {
         public:
             enum class SolutionMethod {
                 SparseLU, Bicgstab
@@ -21,80 +18,72 @@ namespace storm {
                 Ilu, Diagonal, None
             };
             
-            EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method, double precision, uint64_t maximalNumberOfIterations, Preconditioner preconditioner);
-
-            EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            EigenLinearEquationSolverSettings();
             
-            virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            void setSolutionMethod(SolutionMethod const& method);
+            void setPreconditioner(Preconditioner const& preconditioner);
+            void setPrecision(ValueType precision);
+            void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations);
             
-            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            SolutionMethod getSolutionMethod() const;
+            Preconditioner getPreconditioner() const;
+            ValueType getPrecision() const;
+            uint64_t getMaximalNumberOfIterations() const;
             
         private:
-            // A pointer to the original sparse matrix given to this solver.
-            storm::storage::SparseMatrix<ValueType> const* originalA;
-            
-            // The (eigen) matrix associated with this equation solver.
-            std::unique_ptr<Eigen::SparseMatrix<ValueType>> eigenA;
-
-            // The method to use for solving linear equation systems.
             SolutionMethod method;
-            
-            // The required precision for the iterative methods.
-            double precision;
-            
-            // The maximal number of iterations to do before iteration is aborted.
-            uint_fast64_t maximalNumberOfIterations;
-            
-            // The preconditioner to use when solving the linear equation system.
             Preconditioner preconditioner;
+            double precision;
+            uint64_t maximalNumberOfIterations;
         };
         
         template<>
-        class EigenLinearEquationSolver<storm::RationalNumber> : public LinearEquationSolver<storm::RationalNumber> {
+        class EigenLinearEquationSolverSettings<storm::RationalNumber> {
         public:
-            enum class SolutionMethod {
-                SparseLU
-            };
+            EigenLinearEquationSolverSettings();
+        };
+
+        template<>
+        class EigenLinearEquationSolverSettings<storm::RationalFunction> {
+        public:
+            EigenLinearEquationSolverSettings();
+        };
+        
+        /*!
+         * A class that uses the Eigen library to implement the LinearEquationSolver interface.
+         */
+        template<typename ValueType>
+        class EigenLinearEquationSolver : public LinearEquationSolver<ValueType> {
+        public:
+            EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EigenLinearEquationSolverSettings<ValueType> const& settings = EigenLinearEquationSolverSettings<ValueType>());
+            EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EigenLinearEquationSolverSettings<ValueType> const& settings = EigenLinearEquationSolverSettings<ValueType>());
             
-            EigenLinearEquationSolver(storm::storage::SparseMatrix<storm::RationalNumber> const& A, SolutionMethod method = SolutionMethod::SparseLU);
+            virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
-            virtual void solveEquationSystem(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>* multiplyResult = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
-            virtual void performMatrixVectorMultiplication(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const* b, uint_fast64_t n = 1, std::vector<storm::RationalNumber>* multiplyResult = nullptr) const override;
+            EigenLinearEquationSolverSettings<ValueType>& getSettings();
+            EigenLinearEquationSolverSettings<ValueType> const& getSettings() const;
             
         private:
-            // A pointer to the original sparse matrix given to this solver.
-            storm::storage::SparseMatrix<storm::RationalNumber> const* originalA;
-            
             // The (eigen) matrix associated with this equation solver.
-            std::unique_ptr<Eigen::SparseMatrix<storm::RationalNumber>> eigenA;
-            
-            // The method to use for solving linear equation systems.
-            SolutionMethod method;
+            std::unique_ptr<Eigen::SparseMatrix<ValueType>> eigenA;
+
+            // The settings used by the solver.
+            EigenLinearEquationSolverSettings<ValueType> settings;
         };
         
-        template<>
-        class EigenLinearEquationSolver<storm::RationalFunction> : public LinearEquationSolver<storm::RationalFunction> {
+        template<typename ValueType>
+        class EigenLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
         public:
-            enum class SolutionMethod {
-                SparseLU
-            };
-            
-            EigenLinearEquationSolver(storm::storage::SparseMatrix<storm::RationalFunction> const& A, SolutionMethod method = SolutionMethod::SparseLU);
-            
-            virtual void solveEquationSystem(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b, std::vector<storm::RationalFunction>* multiplyResult = nullptr) const override;
-            
-            virtual void performMatrixVectorMultiplication(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const* b, uint_fast64_t n = 1, std::vector<storm::RationalFunction>* multiplyResult = nullptr) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
             
+            EigenLinearEquationSolverSettings<ValueType>& getSettings();
+            EigenLinearEquationSolverSettings<ValueType> const& getSettings() const;
+
         private:
-            // A pointer to the original sparse matrix given to this solver.
-            storm::storage::SparseMatrix<storm::RationalFunction> const* originalA;
-            
-            // The (eigen) matrix associated with this equation solver.
-            std::unique_ptr<Eigen::SparseMatrix<storm::RationalFunction>> eigenA;
-            
-            // The method to use for solving linear equation systems.
-            SolutionMethod method;
+            EigenLinearEquationSolverSettings<ValueType> settings;
         };
     }
 }
\ No newline at end of file
diff --git a/src/solver/EliminationLinearEquationSolver.cpp b/src/solver/EliminationLinearEquationSolver.cpp
index 52a355fbd..e250257ff 100644
--- a/src/solver/EliminationLinearEquationSolver.cpp
+++ b/src/solver/EliminationLinearEquationSolver.cpp
@@ -3,7 +3,6 @@
 #include <numeric>
 
 #include "src/settings/SettingsManager.h"
-#include "src/settings/modules/EliminationSettings.h"
 
 #include "src/solver/stateelimination/StatePriorityQueue.h"
 #include "src/solver/stateelimination/PrioritizedStateEliminator.h"
@@ -20,10 +19,30 @@ namespace storm {
         using namespace storm::utility::stateelimination;
         
         template<typename ValueType>
-        EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) {
-            // Intentionally left empty.
+        EliminationLinearEquationSolverSettings<ValueType>::EliminationLinearEquationSolverSettings() {
+            order = storm::settings::getModule<storm::settings::modules::EliminationSettings>().getEliminationOrder();
         }
         
+        template<typename ValueType>
+        void EliminationLinearEquationSolverSettings<ValueType>::setEliminationOrder(storm::settings::modules::EliminationSettings::EliminationOrder const& order) {
+            this->order = order;
+        }
+        
+        template<typename ValueType>
+        storm::settings::modules::EliminationSettings::EliminationOrder EliminationLinearEquationSolverSettings<ValueType>::getEliminationOrder() const {
+            return order;
+        }
+
+        template<typename ValueType>
+        EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), settings(settings) {
+            // Intentionally left empty.
+        }
+
+        template<typename ValueType>
+        EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), settings(settings) {
+            // Intentionally left empty.
+        }
+
         template<typename ValueType>
         void EliminationLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             STORM_LOG_WARN_COND(multiplyResult == nullptr, "Providing scratch memory will not improve the performance of this solver.");
@@ -36,10 +55,15 @@ namespace storm {
             
             // We need to revert the transformation into an equation system matrix, because the elimination procedure
             // and the distance computation is based on the probability matrix instead.
-            storm::storage::SparseMatrix<ValueType> transitionMatrix(A);
-            
-            transitionMatrix.convertToEquationSystem();
+            storm::storage::SparseMatrix<ValueType> locallyConvertedMatrix;
+            if (localA) {
+                localA->convertToEquationSystem();
+            } else {
+                locallyConvertedMatrix = A;
+                locallyConvertedMatrix.convertToEquationSystem();
+            }
             
+            storm::storage::SparseMatrix<ValueType> const& transitionMatrix = localA ? *localA : locallyConvertedMatrix;
             storm::storage::SparseMatrix<ValueType> backwardTransitions = transitionMatrix.transpose();
 
             // Initialize the solution to the right-hand side of the equation system.
@@ -49,8 +73,8 @@ namespace storm {
             storm::storage::FlexibleSparseMatrix<ValueType> flexibleMatrix(transitionMatrix, false);
             storm::storage::FlexibleSparseMatrix<ValueType> flexibleBackwardTransitions(backwardTransitions, true);
             
-            storm::settings::modules::EliminationSettings::EliminationOrder order = storm::settings::getModule<storm::settings::modules::EliminationSettings>().getEliminationOrder();
             boost::optional<std::vector<uint_fast64_t>> distanceBasedPriorities;
+            auto order = this->getSettings().getEliminationOrder();
             if (eliminationOrderNeedsDistances(order)) {
                 // Since we have no initial states at this point, we determine a representative of every BSCC regarding
                 // the backward transitions, because this means that every row is reachable from this set of rows, which
@@ -69,6 +93,11 @@ namespace storm {
                 auto state = priorityQueue->pop();
                 eliminator.eliminateState(state, false);
             }
+            
+            // After having solved the system, we need to revert the transition system if we kept it local.
+            if (localA) {
+                localA->convertToEquationSystem();
+            };
         }
         
         template<typename ValueType>
@@ -108,10 +137,47 @@ namespace storm {
             }
         }
 
+        template<typename ValueType>
+        EliminationLinearEquationSolverSettings<ValueType>& EliminationLinearEquationSolver<ValueType>::getSettings() {
+            return settings;
+        }
+        
+        template<typename ValueType>
+        EliminationLinearEquationSolverSettings<ValueType> const& EliminationLinearEquationSolver<ValueType>::getSettings() const {
+            return settings;
+        }
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EliminationLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+            return std::make_unique<storm::solver::EliminationLinearEquationSolver<ValueType>>(matrix, settings);
+        }
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EliminationLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
+            return std::make_unique<storm::solver::EliminationLinearEquationSolver<ValueType>>(std::move(matrix), settings);
+        }
+        
+        template<typename ValueType>
+        EliminationLinearEquationSolverSettings<ValueType>& EliminationLinearEquationSolverFactory<ValueType>::getSettings() {
+            return settings;
+        }
+
+        template<typename ValueType>
+        EliminationLinearEquationSolverSettings<ValueType> const& EliminationLinearEquationSolverFactory<ValueType>::getSettings() const {
+            return settings;
+        }
+
+        template class EliminationLinearEquationSolverSettings<double>;
+        template class EliminationLinearEquationSolverSettings<storm::RationalNumber>;
+        template class EliminationLinearEquationSolverSettings<storm::RationalFunction>;
+
         template class EliminationLinearEquationSolver<double>;
         template class EliminationLinearEquationSolver<storm::RationalNumber>;
         template class EliminationLinearEquationSolver<storm::RationalFunction>;
-        
+
+        template class EliminationLinearEquationSolverFactory<double>;
+        template class EliminationLinearEquationSolverFactory<storm::RationalNumber>;
+        template class EliminationLinearEquationSolverFactory<storm::RationalFunction>;
     }
 }
 
diff --git a/src/solver/EliminationLinearEquationSolver.h b/src/solver/EliminationLinearEquationSolver.h
index 25dffa651..dfab51b50 100644
--- a/src/solver/EliminationLinearEquationSolver.h
+++ b/src/solver/EliminationLinearEquationSolver.h
@@ -3,30 +3,65 @@
 
 #include "src/solver/LinearEquationSolver.h"
 
+#include "src/settings/modules/EliminationSettings.h"
+
 namespace storm {
     namespace solver {
+        template<typename ValueType>
+        class EliminationLinearEquationSolverSettings {
+        public:
+            EliminationLinearEquationSolverSettings();
+            
+            void setEliminationOrder(storm::settings::modules::EliminationSettings::EliminationOrder const& order);
+            
+            storm::settings::modules::EliminationSettings::EliminationOrder getEliminationOrder() const;
+
+        private:
+            storm::settings::modules::EliminationSettings::EliminationOrder order;
+        };
+        
         /*!
-         * A class that uses gaussian elimination to implement the LinearEquationSolver interface. In particular
+         * A class that uses gaussian elimination to implement the LinearEquationSolver interface.
          */
         template<typename ValueType>
         class EliminationLinearEquationSolver : public LinearEquationSolver<ValueType> {
         public:
-            
-            /*!
-             * Constructs a linear equation solver.
-             *
-             * @param A The matrix defining the coefficients of the linear equation system.
-             */
-            EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EliminationLinearEquationSolverSettings<ValueType> const& settings = EliminationLinearEquationSolverSettings<ValueType>());
+            EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings = EliminationLinearEquationSolverSettings<ValueType>());
             
             virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
 
+            EliminationLinearEquationSolverSettings<ValueType>& getSettings();
+            EliminationLinearEquationSolverSettings<ValueType> const& getSettings() const;
+            
         private:
-            // A reference to the original matrix used for this equation solver.
+            void initializeSettings();
+            
+            // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
+            // when the solver is destructed.
+            std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
+
+            // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
+            // the reference refers to localA.
             storm::storage::SparseMatrix<ValueType> const& A;
             
+            // The settings used by the solver.
+            EliminationLinearEquationSolverSettings<ValueType> settings;
+        };
+        
+        template<typename ValueType>
+        class EliminationLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
+        public:
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
+            
+            EliminationLinearEquationSolverSettings<ValueType>& getSettings();
+            EliminationLinearEquationSolverSettings<ValueType> const& getSettings() const;
+            
+        private:
+            EliminationLinearEquationSolverSettings<ValueType> settings;
         };
     }
 }
diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp
index 7a0dfdc33..0564938d5 100644
--- a/src/solver/GmmxxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxLinearEquationSolver.cpp
@@ -16,14 +16,8 @@
 
 namespace storm {
     namespace solver {
-        
-        template<typename ValueType>
-        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, bool relative, uint_fast64_t restart) : originalA(&A), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner), relative(relative), restart(restart) {
-            // Intentionally left empty.
-        }
-
         template<typename ValueType>
-        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : originalA(&A), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)) {
+        GmmxxLinearEquationSolverSettings<ValueType>::GmmxxLinearEquationSolverSettings() {
             // Get the settings object to customize linear solving.
             storm::settings::modules::GmmxxEquationSolverSettings const& settings = storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>();
             
@@ -56,40 +50,112 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void GmmxxLinearEquationSolverSettings<ValueType>::setSolutionMethod(SolutionMethod const& method) {
+            this->method = method;
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolverSettings<ValueType>::setPreconditioner(Preconditioner const& preconditioner) {
+            this->preconditioner = preconditioner;
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolverSettings<ValueType>::setPrecision(ValueType precision) {
+            this->precision = precision;
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolverSettings<ValueType>::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) {
+            this->maximalNumberOfIterations = maximalNumberOfIterations;
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolverSettings<ValueType>::setRelativeTerminationCriterion(bool value) {
+            this->relative = value;
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolverSettings<ValueType>::setNumberOfIterationsUntilRestart(uint64_t restart) {
+            this->restart = restart;
+        }
+        
+        template<typename ValueType>
+        typename GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod GmmxxLinearEquationSolverSettings<ValueType>::getSolutionMethod() const {
+            return method;
+        }
+        
+        template<typename ValueType>
+        typename GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner GmmxxLinearEquationSolverSettings<ValueType>::getPreconditioner() const {
+            return preconditioner;
+        }
+        
+        template<typename ValueType>
+        ValueType GmmxxLinearEquationSolverSettings<ValueType>::getPrecision() const {
+            return precision;
+        }
+        
+        template<typename ValueType>
+        uint64_t GmmxxLinearEquationSolverSettings<ValueType>::getMaximalNumberOfIterations() const {
+            return maximalNumberOfIterations;
+        }
+        
+        template<typename ValueType>
+        bool GmmxxLinearEquationSolverSettings<ValueType>::getRelativeTerminationCriterion() const {
+            return relative;
+        }
+        
+        template<typename ValueType>
+        uint64_t GmmxxLinearEquationSolverSettings<ValueType>::getNumberOfIterationsUntilRestart() const {
+            return restart;
+        }
+        
+        template<typename ValueType>
+        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), settings(settings) {
+            // Intentionally left empty.
+        }
+
+        template<typename ValueType>
+        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(*localA)), settings(settings) {
+            // Intentionally left empty.
+        }
+        
         template<typename ValueType>
         void GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
-            STORM_LOG_INFO("Using method '" << methodToString() << "' with preconditioner '" << preconditionerToString() << "' (max. " << maximalNumberOfIterations << " iterations).");
-            if (method == SolutionMethod::Jacobi && preconditioner != Preconditioner::None) {
+            auto method = this->getSettings().getSolutionMethod();
+            auto preconditioner = this->getSettings().getPreconditioner();
+            STORM_LOG_INFO("Using method '" << method << "' with preconditioner '" << preconditioner << "' (max. " << this->getSettings().getMaximalNumberOfIterations() << " iterations).");
+            if (method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi && preconditioner != GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::None) {
                 STORM_LOG_WARN("Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored.");
             }
             
-            if (method == SolutionMethod::Bicgstab || method == SolutionMethod::Qmr || method == SolutionMethod::Gmres) {
+            if (method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Bicgstab || method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Qmr || method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Gmres) {
                 // Prepare an iteration object that determines the accuracy and the maximum number of iterations.
-                gmm::iteration iter(precision, 0, maximalNumberOfIterations);
+                gmm::iteration iter(this->getSettings().getPrecision(), 0, this->getSettings().getMaximalNumberOfIterations());
                 
-                if (method == SolutionMethod::Bicgstab) {
-                    if (preconditioner == Preconditioner::Ilu) {
+                if (method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Bicgstab) {
+                    if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::Ilu) {
                         gmm::bicgstab(*gmmxxMatrix, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), iter);
-                    } else if (preconditioner == Preconditioner::Diagonal) {
+                    } else if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::Diagonal) {
                         gmm::bicgstab(*gmmxxMatrix, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), iter);
-                    } else if (preconditioner == Preconditioner::None) {
+                    } else if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::None) {
                         gmm::bicgstab(*gmmxxMatrix, x, b, gmm::identity_matrix(), iter);
                     }
-                } else if (method == SolutionMethod::Qmr) {
-                    if (preconditioner == Preconditioner::Ilu) {
+                } else if (method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Qmr) {
+                    if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::Ilu) {
                         gmm::qmr(*gmmxxMatrix, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), iter);
-                    } else if (preconditioner == Preconditioner::Diagonal) {
+                    } else if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::Diagonal) {
                         gmm::qmr(*gmmxxMatrix, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), iter);
-                    } else if (preconditioner == Preconditioner::None) {
+                    } else if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::None) {
                         gmm::qmr(*gmmxxMatrix, x, b, gmm::identity_matrix(), iter);
                     }
-                } else if (method == SolutionMethod::Gmres) {
-                    if (preconditioner == Preconditioner::Ilu) {
-                        gmm::gmres(*gmmxxMatrix, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), restart, iter);
-                    } else if (preconditioner == Preconditioner::Diagonal) {
-                        gmm::gmres(*gmmxxMatrix, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), restart, iter);
-                    } else if (preconditioner == Preconditioner::None) {
-                        gmm::gmres(*gmmxxMatrix, x, b, gmm::identity_matrix(), restart, iter);
+                } else if (method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Gmres) {
+                    if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::Ilu) {
+                        gmm::gmres(*gmmxxMatrix, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), this->getSettings().getNumberOfIterationsUntilRestart(), iter);
+                    } else if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::Diagonal) {
+                        gmm::gmres(*gmmxxMatrix, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), this->getSettings().getNumberOfIterationsUntilRestart(), iter);
+                    } else if (preconditioner == GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::None) {
+                        gmm::gmres(*gmmxxMatrix, x, b, gmm::identity_matrix(), this->getSettings().getNumberOfIterationsUntilRestart(), iter);
                     }
                 }
                 
@@ -99,11 +165,11 @@ namespace storm {
                 } else {
                     STORM_LOG_WARN("Iterative solver did not converge.");
                 }
-            } else if (method == SolutionMethod::Jacobi) {
-                uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(*originalA, x, b, multiplyResult);
+            } else if (method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
+                uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b, multiplyResult);
                 
                 // Check if the solver converged and issue a warning otherwise.
-                if (iterations < maximalNumberOfIterations) {
+                if (iterations < this->getSettings().getMaximalNumberOfIterations()) {
                     STORM_LOG_INFO("Iterative solver converged after " << iterations << " iterations.");
                 } else {
                     STORM_LOG_WARN("Iterative solver did not converge.");
@@ -173,14 +239,14 @@ namespace storm {
             uint_fast64_t iterationCount = 0;
             bool converged = false;
             
-            while (!converged && iterationCount < maximalNumberOfIterations && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
+            while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
                 // Compute D^-1 * (b - LU * x) and store result in nextX.
                 gmm::mult(*gmmLU, *currentX, tmpX);
                 gmm::add(b, gmm::scaled(tmpX, -storm::utility::one<ValueType>()), tmpX);
                 storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
                 
                 // Now check if the process already converged within our precision.
-                converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, precision, relative);
+                converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion());
 
                 // Swap the two pointers as a preparation for the next iteration.
                 std::swap(nextX, currentX);
@@ -204,27 +270,39 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::string GmmxxLinearEquationSolver<ValueType>::methodToString() const {
-            switch (method) {
-                case SolutionMethod::Bicgstab: return "bicgstab";
-                case SolutionMethod::Qmr: return "qmr";
-                case SolutionMethod::Gmres: return "gmres";
-                case SolutionMethod::Jacobi: return "jacobi";
-                default: return "invalid";
-            }
+        GmmxxLinearEquationSolverSettings<ValueType>& GmmxxLinearEquationSolver<ValueType>::getSettings() {
+            return settings;
+        }
+
+        template<typename ValueType>
+        GmmxxLinearEquationSolverSettings<ValueType> const& GmmxxLinearEquationSolver<ValueType>::getSettings() const {
+            return settings;
         }
         
         template<typename ValueType>
-        std::string GmmxxLinearEquationSolver<ValueType>::preconditionerToString() const {
-            switch (preconditioner) {
-                case Preconditioner::Ilu: return "ilu";
-                case Preconditioner::Diagonal: return "diagonal";
-                case Preconditioner::None: return "none";
-                default: return "invalid";
-            }
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+            return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(matrix, settings);
+        }
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
+            return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(std::move(matrix), settings);
+        }
+        
+        template<typename ValueType>
+        GmmxxLinearEquationSolverSettings<ValueType>& GmmxxLinearEquationSolverFactory<ValueType>::getSettings() {
+            return settings;
+        }
+        
+        template<typename ValueType>
+        GmmxxLinearEquationSolverSettings<ValueType> const& GmmxxLinearEquationSolverFactory<ValueType>::getSettings() const {
+            return settings;
         }
         
         // Explicitly instantiate the solver.
+        template class GmmxxLinearEquationSolverSettings<double>;
         template class GmmxxLinearEquationSolver<double>;
+        template class GmmxxLinearEquationSolverFactory<double>;
+        
     }
 }
\ No newline at end of file
diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h
index a6c402f89..89448e4f3 100644
--- a/src/solver/GmmxxLinearEquationSolver.h
+++ b/src/solver/GmmxxLinearEquationSolver.h
@@ -1,6 +1,8 @@
 #ifndef STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_
 #define STORM_SOLVER_GMMXXLINEAREQUATIONSOLVER_H_
 
+#include <ostream>
+
 #include "src/utility/gmm.h"
 
 #include "LinearEquationSolver.h"
@@ -8,48 +10,92 @@
 namespace storm {
     namespace solver {
 
-        /*!
-         * A class that uses the gmm++ library to implement the LinearEquationSolver interface.
-         */
         template<typename ValueType>
-        class GmmxxLinearEquationSolver : public LinearEquationSolver<ValueType> {
+        class GmmxxLinearEquationSolverSettings {
         public:
             // An enumeration specifying the available preconditioners.
             enum class Preconditioner {
                 Ilu, Diagonal, None
             };
+
             
+            friend std::ostream& operator<<(std::ostream& out, Preconditioner const& preconditioner) {
+                switch (preconditioner) {
+                    case GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::Ilu: out << "ilu"; break;
+                    case GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::Diagonal: out << "diagonal"; break;
+                    case GmmxxLinearEquationSolverSettings<ValueType>::Preconditioner::None: out << "none"; break;
+                }
+                return out;
+            }
+
             // An enumeration specifying the available solution methods.
             enum class SolutionMethod {
                 Bicgstab, Qmr, Gmres, Jacobi
             };
             
-            /*!
-             * Constructs a linear equation solver with the given parameters.
-             *
-             * @param A The matrix defining the coefficients of the linear equation system.
-             * @param method The method to use for linear equation solving.
-             * @param precision The precision to use for convergence detection.
-             * @param maximalNumberOfIterations The maximal number of iterations do perform before iteration is aborted.
-             * @param preconditioner The preconditioner to use when solving linear equation systems.
-             * @param relative If set, the relative error rather than the absolute error is considered for convergence
-             * detection.
-             * @param restart An optional argument that specifies after how many iterations restarted methods are
-             * supposed to actually to a restart.
-             */
-            GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, bool relative = true, uint_fast64_t restart = 0);
+            friend std::ostream& operator<<(std::ostream& out, SolutionMethod const& method) {
+                switch (method) {
+                    case GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Bicgstab: out << "BiCGSTAB"; break;
+                    case GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Qmr: out << "QMR"; break;
+                    case GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Gmres: out << "GMRES"; break;
+                    case GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi: out << "Jacobi"; break;
+                }
+                return out;
+            }
+
+            GmmxxLinearEquationSolverSettings();
+            
+            void setSolutionMethod(SolutionMethod const& method);
+            void setPreconditioner(Preconditioner const& preconditioner);
+            void setPrecision(ValueType precision);
+            void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations);
+            void setRelativeTerminationCriterion(bool value);
+            void setNumberOfIterationsUntilRestart(uint64_t restart);
+         
+            SolutionMethod getSolutionMethod() const;
+            Preconditioner getPreconditioner() const;
+            ValueType getPrecision() const;
+            uint64_t getMaximalNumberOfIterations() const;
+            bool getRelativeTerminationCriterion() const;
+            uint64_t getNumberOfIterationsUntilRestart() const;
+            
+        private:
+            // The method to use for solving linear equation systems.
+            SolutionMethod method;
+            
+            // The required precision for the iterative methods.
+            double precision;
+            
+            // The maximal number of iterations to do before iteration is aborted.
+            uint_fast64_t maximalNumberOfIterations;
+            
+            // The preconditioner to use when solving the linear equation system.
+            Preconditioner preconditioner;
+            
+            // Sets whether the relative or absolute error is to be considered for convergence detection. Note that this
+            // only applies to the Jacobi method for this solver.
+            bool relative;
+            
+            // A restart value that determines when restarted methods shall do so.
+            uint_fast64_t restart;
+        };
+        
+        /*!
+         * A class that uses the gmm++ library to implement the LinearEquationSolver interface.
+         */
+        template<typename ValueType>
+        class GmmxxLinearEquationSolver : public LinearEquationSolver<ValueType> {
+        public:
+            GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings = GmmxxLinearEquationSolverSettings<ValueType>());
+            GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings = GmmxxLinearEquationSolverSettings<ValueType>());
             
-            /*!
-             * Constructs a linear equation solver with parameters being set according to the settings object.
-             *
-             * @param A The matrix defining the coefficients of the linear equation system.
-             */
-            GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
-                        
             virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
+            GmmxxLinearEquationSolverSettings<ValueType>& getSettings();
+            GmmxxLinearEquationSolverSettings<ValueType> const& getSettings() const;
+            
         private:
             /*!
              * Solves the linear equation system A*x = b given by the parameters using the Jacobi method.
@@ -65,44 +111,32 @@ namespace storm {
              */
             uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const;
 
-            /*!
-             * Retrieves the string representation of the solution method associated with this solver.
-             *
-             * @return The string representation of the solution method associated with this solver.
-             */
-            std::string methodToString() const;
-            
-            /*!
-             * Retrieves the string representation of the preconditioner associated with this solver.
-             *
-             * @return The string representation of the preconditioner associated with this solver.
-             */
-            std::string preconditionerToString() const;
-            
-            // A pointer to the original sparse matrix given to this solver.
-            storm::storage::SparseMatrix<ValueType> const* originalA;
+            // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
+            // when the solver is destructed.
+            std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
+
+            // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
+            // the reference refers to localA.
+            storm::storage::SparseMatrix<ValueType> const& A;
             
             // The (gmm++) matrix associated with this equation solver.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix;
             
-            // The method to use for solving linear equation systems.
-            SolutionMethod method;
-
-            // The required precision for the iterative methods.
-            double precision;
-
-            // The maximal number of iterations to do before iteration is aborted.
-            uint_fast64_t maximalNumberOfIterations;
-
-            // The preconditioner to use when solving the linear equation system.
-            Preconditioner preconditioner;
-
-            // Sets whether the relative or absolute error is to be considered for convergence detection. Note that this
-            // only applies to the Jacobi method for this solver.
-            bool relative;
-
-            // A restart value that determines when restarted methods shall do so.
-            uint_fast64_t restart;
+            // The settings used by the solver.
+            GmmxxLinearEquationSolverSettings<ValueType> settings;
+        };
+        
+        template<typename ValueType>
+        class GmmxxLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
+        public:
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
+            
+            GmmxxLinearEquationSolverSettings<ValueType>& getSettings();
+            GmmxxLinearEquationSolverSettings<ValueType> const& getSettings() const;
+            
+        private:
+            GmmxxLinearEquationSolverSettings<ValueType> settings;
         };
 
     } // namespace solver
diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
index e616f8d62..02835670f 100644
--- a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
@@ -122,14 +122,15 @@ namespace storm {
 					storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(scheduler, true);
 					submatrix.convertToEquationSystem();
 					
-					GmmxxLinearEquationSolver<ValueType> gmmLinearEquationSolver(submatrix, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, this->precision, this->maximalNumberOfIterations, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Ilu, this->relative, 50);
+					GmmxxLinearEquationSolver<ValueType> gmmxxLinearEquationSolver(submatrix);
+                    
 					storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, rowGroupIndices, b);
 
 					// Copy X since we will overwrite it
 					std::copy(currentX->begin(), currentX->end(), newX->begin());
 
 					// Solve the resulting linear equation system
-					gmmLinearEquationSolver.solveEquationSystem(*newX, subB, &deterministicMultiplyResult);
+					gmmxxLinearEquationSolver.solveEquationSystem(*newX, subB, &deterministicMultiplyResult);
 
 					// Compute x' = A*x + b. This step is necessary to allow the choosing of the optimal policy for the next iteration.
 					gmm::mult(*gmmxxMatrix, *newX, *multiplyResult);
diff --git a/src/solver/LinearEquationSolver.cpp b/src/solver/LinearEquationSolver.cpp
new file mode 100644
index 000000000..801855765
--- /dev/null
+++ b/src/solver/LinearEquationSolver.cpp
@@ -0,0 +1,87 @@
+#include "src/solver/LinearEquationSolver.h"
+
+#include "src/solver/SolverSelectionOptions.h"
+
+#include "src/solver/GmmxxLinearEquationSolver.h"
+#include "src/solver/NativeLinearEquationSolver.h"
+#include "src/solver/EigenLinearEquationSolver.h"
+#include "src/solver/EliminationLinearEquationSolver.h"
+
+#include "src/settings/SettingsManager.h"
+#include "src/settings/modules/MarkovChainSettings.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
+            return create(matrix);
+        }
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GeneralLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+            return selectSolver(matrix);
+        }
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GeneralLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
+            return selectSolver(std::move(matrix));
+        }
+        
+        template<typename ValueType>
+        template<typename MatrixType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GeneralLinearEquationSolverFactory<ValueType>::selectSolver(MatrixType&& matrix) const {
+            storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::MarkovChainSettings>().getEquationSolver();
+            switch (equationSolver) {
+                case storm::solver::EquationSolverType::Gmmxx: return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
+                case storm::solver::EquationSolverType::Native: return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
+                case storm::solver::EquationSolverType::Eigen: return std::make_unique<storm::solver::EigenLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
+                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::EliminationLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
+                default: return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
+            }
+        }
+        
+        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> GeneralLinearEquationSolverFactory<storm::RationalNumber>::create(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix) const {
+            return selectSolver(matrix);
+        }
+        
+        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> GeneralLinearEquationSolverFactory<storm::RationalNumber>::create(storm::storage::SparseMatrix<storm::RationalNumber>&& matrix) const {
+            return selectSolver(std::move(matrix));
+        }
+        
+        template<typename MatrixType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> GeneralLinearEquationSolverFactory<storm::RationalNumber>::selectSolver(MatrixType&& matrix) const {
+            storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::MarkovChainSettings>().getEquationSolver();
+            switch (equationSolver) {
+                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::EliminationLinearEquationSolver<storm::RationalNumber>>(matrix);
+                default: return std::make_unique<storm::solver::EigenLinearEquationSolver<storm::RationalNumber>>(matrix);
+            }
+        }
+        
+        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> GeneralLinearEquationSolverFactory<storm::RationalFunction>::create(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix) const {
+            return selectSolver(matrix);
+        }
+        
+        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> GeneralLinearEquationSolverFactory<storm::RationalFunction>::create(storm::storage::SparseMatrix<storm::RationalFunction>&& matrix) const {
+            return selectSolver(std::move(matrix));
+        }
+        
+        template<typename MatrixType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> GeneralLinearEquationSolverFactory<storm::RationalFunction>::selectSolver(MatrixType&& matrix) const {
+            storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::MarkovChainSettings>().getEquationSolver();
+            switch (equationSolver) {
+                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::EliminationLinearEquationSolver<storm::RationalFunction>>(matrix);
+                default: return std::make_unique<storm::solver::EigenLinearEquationSolver<storm::RationalFunction>>(matrix);
+            }
+        }
+        
+        template class LinearEquationSolverFactory<double>;
+        template class LinearEquationSolverFactory<storm::RationalNumber>;
+        template class LinearEquationSolverFactory<storm::RationalFunction>;
+        
+        template class GeneralLinearEquationSolverFactory<double>;
+        template class GeneralLinearEquationSolverFactory<storm::RationalNumber>;
+        template class GeneralLinearEquationSolverFactory<storm::RationalFunction>;
+
+    }
+}
\ No newline at end of file
diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h
index e85fe3cce..833a48229 100644
--- a/src/solver/LinearEquationSolver.h
+++ b/src/solver/LinearEquationSolver.h
@@ -48,6 +48,60 @@ namespace storm {
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
         };
         
+        template<typename ValueType>
+        class LinearEquationSolverFactory {
+        public:
+            /*!
+             * Creates a new linear equation solver instance with the given matrix.
+             *
+             * @param matrix The matrix that defines the equation system.
+             * @return A pointer to the newly created solver.
+             */
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const = 0;
+            
+            /*!
+             * Creates a new linear equation solver instance with the given matrix. The caller gives up posession of the
+             * matrix by calling this function.
+             *
+             * @param matrix The matrix that defines the equation system.
+             * @return A pointer to the newly created solver.
+             */
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const;
+        };
+        
+        template<typename ValueType>
+        class GeneralLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
+        public:
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
+            
+        private:
+            template<typename MatrixType>
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> selectSolver(MatrixType&& matrix) const;
+        };
+        
+        template<>
+        class GeneralLinearEquationSolverFactory<storm::RationalNumber> : public LinearEquationSolverFactory<storm::RationalNumber> {
+        public:
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> create(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> create(storm::storage::SparseMatrix<storm::RationalNumber>&& matrix) const override;
+            
+        private:
+            template<typename MatrixType>
+            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> selectSolver(MatrixType&& matrix) const;
+        };
+        
+        template<>
+        class GeneralLinearEquationSolverFactory<storm::RationalFunction> : public LinearEquationSolverFactory<storm::RationalFunction> {
+        public:
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> create(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> create(storm::storage::SparseMatrix<storm::RationalFunction>&& matrix) const override;
+            
+        private:
+            template<typename MatrixType>
+            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> selectSolver(MatrixType&& matrix) const;
+        };
+        
     } // namespace solver
 } // namespace storm
 
diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp
index 9c402255b..bcfaf4782 100644
--- a/src/solver/NativeLinearEquationSolver.cpp
+++ b/src/solver/NativeLinearEquationSolver.cpp
@@ -7,31 +7,97 @@
 
 #include "src/utility/vector.h"
 #include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/InvalidSettingsException.h"
 
 namespace storm {
     namespace solver {
-        
-        template<typename ValueType>
-        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), method(method), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
-            // Intentionally left empty.
-        }
-        
+
         template<typename ValueType>
-        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSolutionMethod method) : A(A), method(method) {
-            // Get the settings object to customize linear solving.
+        NativeLinearEquationSolverSettings<ValueType>::NativeLinearEquationSolverSettings() {
             storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>();
             
-            // Get appropriate settings.
+            storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod methodAsSetting = settings.getLinearEquationSystemMethod();
+            if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel) {
+                method = SolutionMethod::GaussSeidel;
+            } else if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Jacobi) {
+                method = SolutionMethod::Jacobi;
+            } else if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::SOR) {
+                method = SolutionMethod::SOR;
+            } else {
+                STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "The selected solution technique is invalid for this solver.");
+            }
+            
             maximalNumberOfIterations = settings.getMaximalIterationCount();
             precision = settings.getPrecision();
             relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
+            omega = storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getOmega();
         }
-                
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolverSettings<ValueType>::setSolutionMethod(SolutionMethod const& method) {
+            this->method = method;
+        }
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolverSettings<ValueType>::setPrecision(ValueType precision) {
+            this->precision = precision;
+        }
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolverSettings<ValueType>::setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations) {
+            this->maximalNumberOfIterations = maximalNumberOfIterations;
+        }
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolverSettings<ValueType>::setRelativeTerminationCriterion(bool value) {
+            this->relative = value;
+        }
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolverSettings<ValueType>::setOmega(ValueType omega) {
+            this->omega = omega;
+        }
+        
+        template<typename ValueType>
+        typename NativeLinearEquationSolverSettings<ValueType>::SolutionMethod NativeLinearEquationSolverSettings<ValueType>::getSolutionMethod() const {
+            return method;
+        }
+        
+        template<typename ValueType>
+        ValueType NativeLinearEquationSolverSettings<ValueType>::getPrecision() const {
+            return precision;
+        }
+        
+        template<typename ValueType>
+        uint64_t NativeLinearEquationSolverSettings<ValueType>::getMaximalNumberOfIterations() const {
+            return maximalNumberOfIterations;
+        }
+        
+        template<typename ValueType>
+        uint64_t NativeLinearEquationSolverSettings<ValueType>::getRelativeTerminationCriterion() const {
+            return relative;
+        }
+        
+        template<typename ValueType>
+        ValueType NativeLinearEquationSolverSettings<ValueType>::getOmega() const {
+            return omega;
+        }
+        
+        template<typename ValueType>
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), settings(settings) {
+            // Intentionally left empty.
+        }
+
+        template<typename ValueType>
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), settings(settings) {
+            // Intentionally left empty.
+        }
+        
         template<typename ValueType>
         void NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
-            if (method == NativeLinearEquationSolverSolutionMethod::SOR || method == NativeLinearEquationSolverSolutionMethod::GaussSeidel) {
+            if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::GaussSeidel) {
                 // Define the omega used for SOR.
-                ValueType omega = method == NativeLinearEquationSolverSolutionMethod::SOR ? storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getOmega() : storm::utility::one<ValueType>();
+                ValueType omega = this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::SOR ? this->getSettings().getOmega() : storm::utility::one<ValueType>();
                 
                 // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
                 bool tmpXProvided = true;
@@ -47,11 +113,11 @@ namespace storm {
                 uint_fast64_t iterationCount = 0;
                 bool converged = false;
                 
-                while (!converged && iterationCount < maximalNumberOfIterations) {
+                while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations()) {
                     A.performSuccessiveOverRelaxationStep(omega, x, b);
                     
                     // Now check if the process already converged within our precision.
-                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(x, *tmpX, static_cast<ValueType>(precision), relative) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x));
+                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(x, *tmpX, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x));
                     
                     // If we did not yet converge, we need to copy the contents of x to *tmpX.
                     if (!converged) {
@@ -87,7 +153,7 @@ namespace storm {
                 uint_fast64_t iterationCount = 0;
                 bool converged = false;
                 
-                while (!converged && iterationCount < maximalNumberOfIterations && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
+                while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
                     // Compute D^-1 * (b - LU * x) and store result in nextX.
                     jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
                     storm::utility::vector::subtractVectors(b, tmpX, tmpX);
@@ -97,7 +163,7 @@ namespace storm {
                     std::swap(nextX, currentX);
                     
                     // Now check if the process already converged within our precision.
-                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(precision), relative);
+                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion());
                     
                     // Increase iteration count so we can abort if convergence is too slow.
                     ++iterationCount;
@@ -152,20 +218,40 @@ namespace storm {
                 delete copyX;
             }
         }
-
-
+        
         template<typename ValueType>
-        std::string NativeLinearEquationSolver<ValueType>::methodToString() const {
-            switch (method) {
-                case NativeLinearEquationSolverSolutionMethod::Jacobi: return "jacobi";
-                case NativeLinearEquationSolverSolutionMethod::GaussSeidel: return "gauss-seidel";
-                case NativeLinearEquationSolverSolutionMethod::SOR: return "sor";
-                default: return "invalid";
-            }
+        NativeLinearEquationSolverSettings<ValueType>& NativeLinearEquationSolver<ValueType>::getSettings() {
+            return settings;
+        }
+        
+        template<typename ValueType>
+        NativeLinearEquationSolverSettings<ValueType> const& NativeLinearEquationSolver<ValueType>::getSettings() const {
+            return settings;
+        }
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+            return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(matrix, settings);
+        }
+        
+        template<typename ValueType>
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
+            return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(std::move(matrix), settings);
+        }
+        
+        template<typename ValueType>
+        NativeLinearEquationSolverSettings<ValueType>& NativeLinearEquationSolverFactory<ValueType>::getSettings() {
+            return settings;
+        }
+        
+        template<typename ValueType>
+        NativeLinearEquationSolverSettings<ValueType> const& NativeLinearEquationSolverFactory<ValueType>::getSettings() const {
+            return settings;
         }
         
         // Explicitly instantiate the linear equation solver.
+        template class NativeLinearEquationSolverSettings<double>;
         template class NativeLinearEquationSolver<double>;
-		template class NativeLinearEquationSolver<float>;
+        template class NativeLinearEquationSolverFactory<double>;
     }
 }
\ No newline at end of file
diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h
index 9dabd4e65..98abc4264 100644
--- a/src/solver/NativeLinearEquationSolver.h
+++ b/src/solver/NativeLinearEquationSolver.h
@@ -8,66 +8,75 @@
 namespace storm {
     namespace solver {
         
-        // An enumeration specifying the available solution methods.
-        enum class NativeLinearEquationSolverSolutionMethod {
-            Jacobi, GaussSeidel, SOR
+        template<typename ValueType>
+        class NativeLinearEquationSolverSettings {
+        public:
+            enum class SolutionMethod {
+                Jacobi, GaussSeidel, SOR
+            };
+
+            NativeLinearEquationSolverSettings();
+            
+            void setSolutionMethod(SolutionMethod const& method);
+            void setPrecision(ValueType precision);
+            void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations);
+            void setRelativeTerminationCriterion(bool value);
+            void setOmega(ValueType omega);
+            
+            SolutionMethod getSolutionMethod() const;
+            ValueType getPrecision() const;
+            uint64_t getMaximalNumberOfIterations() const;
+            uint64_t getRelativeTerminationCriterion() const;
+            ValueType getOmega() const;
+            
+        private:
+            SolutionMethod method;
+            double precision;
+            bool relative;
+            uint_fast64_t maximalNumberOfIterations;
+            ValueType omega;
         };
         
-        std::ostream& operator<<(std::ostream& out, NativeLinearEquationSolverSolutionMethod const& method);
-        
         /*!
          * A class that uses Storm's native matrix operations to implement the LinearEquationSolver interface.
          */
         template<typename ValueType>
         class NativeLinearEquationSolver : public LinearEquationSolver<ValueType> {
         public:
+            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings = NativeLinearEquationSolverSettings<ValueType>());
+            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings = NativeLinearEquationSolverSettings<ValueType>());
             
-            /*!
-             * Constructs a linear equation solver with parameters being set according to the settings object.
-             *
-             * @param A The matrix defining the coefficients of the linear equation system.
-             * @param method The method to use for solving linear equations.
-             */
-            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSolutionMethod method = NativeLinearEquationSolverSolutionMethod::GaussSeidel);
-
-            /*!
-             * Constructs a linear equation solver with the given parameters.
-             *
-             * @param A The matrix defining the coefficients of the linear equation system.
-             * @param method The method to use for linear equation solving.
-             * @param precision The precision to use for convergence detection.
-             * @param maximalNumberOfIterations The maximal number of iterations do perform before iteration is aborted.
-             * @param relative If set, the relative error rather than the absolute error is considered for convergence
-             * detection.
-             */
-            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
-                        
             virtual void solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
             virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
 
+            NativeLinearEquationSolverSettings<ValueType>& getSettings();
+            NativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
+
         private:
-            /*!
-             * Retrieves the string representation of the solution method associated with this solver.
-             *
-             * @return The string representation of the solution method associated with this solver.
-             */
-            std::string methodToString() const;
+            // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
+            // when the solver is destructed.
+            std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
             
-            // A reference to the matrix the gives the coefficients of the linear equation system.
+            // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
+            // the reference refers to localA.
             storm::storage::SparseMatrix<ValueType> const& A;
+                        
+            // The settings used by the solver.
+            NativeLinearEquationSolverSettings<ValueType> settings;
+        };
+        
+        template<typename ValueType>
+        class NativeLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
+        public:
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
             
-            // The method to use for solving linear equation systems.
-            NativeLinearEquationSolverSolutionMethod method;
-            
-            // The required precision for the iterative methods.
-            double precision;
-            
-            // Sets whether the relative or absolute error is to be considered for convergence detection.
-            bool relative;
+            NativeLinearEquationSolverSettings<ValueType>& getSettings();
+            NativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
             
-            // The maximal number of iterations to do before iteration is aborted.
-            uint_fast64_t maximalNumberOfIterations;
+        private:
+            NativeLinearEquationSolverSettings<ValueType> settings;
         };
     }
 }
diff --git a/src/solver/NativeMinMaxLinearEquationSolver.cpp b/src/solver/NativeMinMaxLinearEquationSolver.cpp
index 544840898..26d34f991 100644
--- a/src/solver/NativeMinMaxLinearEquationSolver.cpp
+++ b/src/solver/NativeMinMaxLinearEquationSolver.cpp
@@ -205,6 +205,5 @@ namespace storm {
         
         // Explicitly instantiate the solver.
         template class NativeMinMaxLinearEquationSolver<double>;
-		template class NativeMinMaxLinearEquationSolver<float>;
     } // namespace solver
 } // namespace storm
diff --git a/src/solver/SolveGoal.cpp b/src/solver/SolveGoal.cpp
index 642077c05..24ca63ea8 100644
--- a/src/solver/SolveGoal.cpp
+++ b/src/solver/SolveGoal.cpp
@@ -32,14 +32,14 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(BoundedGoal<ValueType> const& goal, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix) {
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(BoundedGoal<ValueType> const& goal, storm::solver::LinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix) {
             std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = factory.create(matrix);
             solver->setTerminationCondition(std::make_unique<TerminateIfFilteredExtremumExceedsThreshold<double>>(goal.relevantValues(), goal.thresholdValue(), goal.boundIsStrict(), goal.minimize()));
             return solver;
         }
         
         template<typename ValueType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(SolveGoal const& goal, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix) {
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(SolveGoal const& goal, storm::solver::LinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix) {
             if (goal.isBounded()) {
                 return configureLinearEquationSolver(static_cast<BoundedGoal<ValueType> const&>(goal), factory, matrix);
             }
@@ -48,8 +48,8 @@ namespace storm {
     
         template std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<double>> configureMinMaxLinearEquationSolver(BoundedGoal<double> const& goal, storm::utility::solver::MinMaxLinearEquationSolverFactory<double> const& factory, storm::storage::SparseMatrix<double> const& matrix);
         template std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<double>> configureMinMaxLinearEquationSolver(SolveGoal const& goal, storm::utility::solver::MinMaxLinearEquationSolverFactory<double> const& factory, storm::storage::SparseMatrix<double> const&  matrix);
-        template std::unique_ptr<storm::solver::LinearEquationSolver<double>> configureLinearEquationSolver(BoundedGoal<double> const& goal, storm::utility::solver::LinearEquationSolverFactory<double> const& factory, storm::storage::SparseMatrix<double> const&  matrix);
-        template std::unique_ptr<storm::solver::LinearEquationSolver<double>> configureLinearEquationSolver(SolveGoal const& goal, storm::utility::solver::LinearEquationSolverFactory<double> const& factory, storm::storage::SparseMatrix<double> const&  matrix);
+        template std::unique_ptr<storm::solver::LinearEquationSolver<double>> configureLinearEquationSolver(BoundedGoal<double> const& goal, storm::solver::LinearEquationSolverFactory<double> const& factory, storm::storage::SparseMatrix<double> const&  matrix);
+        template std::unique_ptr<storm::solver::LinearEquationSolver<double>> configureLinearEquationSolver(SolveGoal const& goal, storm::solver::LinearEquationSolverFactory<double> const& factory, storm::storage::SparseMatrix<double> const&  matrix);
         
     }
 }
diff --git a/src/solver/SolveGoal.h b/src/solver/SolveGoal.h
index fb28dfbb3..2dab69c4f 100644
--- a/src/solver/SolveGoal.h
+++ b/src/solver/SolveGoal.h
@@ -8,6 +8,8 @@
 #include "src/logic/Bound.h"
 #include "src/storage/BitVector.h"
 
+#include "src/solver/LinearEquationSolver.h"
+
 namespace storm {
     namespace storage {
         template<typename ValueType> class SparseMatrix;
@@ -101,10 +103,10 @@ namespace storm {
         std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> configureMinMaxLinearEquationSolver(SolveGoal const& goal, storm::utility::solver::MinMaxLinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix);
 
         template<typename ValueType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(BoundedGoal<ValueType> const& goal, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix);
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(BoundedGoal<ValueType> const& goal, storm::solver::LinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix);
         
         template<typename ValueType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(SolveGoal const& goal, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix);
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(SolveGoal const& goal, storm::solver::LinearEquationSolverFactory<ValueType> const& factory, storm::storage::SparseMatrix<ValueType> const& matrix);
 
     }
 }
diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
index 0140fec1d..2f6e51de1 100644
--- a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
+++ b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
@@ -429,6 +429,5 @@ namespace storm {
 
         // Explicitly instantiate the solver.
 		template class TopologicalMinMaxLinearEquationSolver<double>;
-		template class TopologicalMinMaxLinearEquationSolver<float>;
     } // namespace solver
 } // namespace storm
diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp
index 242078c9e..cafe63565 100644
--- a/src/utility/solver.cpp
+++ b/src/utility/solver.cpp
@@ -45,119 +45,7 @@ namespace storm {
             std::unique_ptr<storm::solver::SymbolicGameSolver<Type, ValueType>> SymbolicGameSolverFactory<Type, ValueType>::create(storm::dd::Add<Type, ValueType> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) const {
                 return std::unique_ptr<storm::solver::SymbolicGameSolver<Type, ValueType>>(new storm::solver::SymbolicGameSolver<Type, ValueType>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, player1Variables, player2Variables));
             }
-            
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
-                return create(matrix);
-            }
-            
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GeneralLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
-                return selectSolver(matrix);
-            }
-            
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GeneralLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
-                return selectSolver(std::move(matrix));
-            }
-            
-            template<typename ValueType>
-            template<typename MatrixType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GeneralLinearEquationSolverFactory<ValueType>::selectSolver(MatrixType&& matrix) const {
-                storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::MarkovChainSettings>().getEquationSolver();
-                switch (equationSolver) {
-                    case storm::solver::EquationSolverType::Gmmxx: return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
-                    case storm::solver::EquationSolverType::Native: return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
-                    case storm::solver::EquationSolverType::Eigen: return std::make_unique<storm::solver::EigenLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
-                    case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::EliminationLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
-                    default: return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(std::forward<MatrixType>(matrix));
-                }
-            }
-
-            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> GeneralLinearEquationSolverFactory<storm::RationalNumber>::create(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix) const {
-                return selectSolver(matrix);
-            }
 
-            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> GeneralLinearEquationSolverFactory<storm::RationalNumber>::create(storm::storage::SparseMatrix<storm::RationalNumber>&& matrix) const {
-                return selectSolver(std::move(matrix));
-            }
-            
-            template<typename MatrixType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> GeneralLinearEquationSolverFactory<storm::RationalNumber>::selectSolver(MatrixType&& matrix) const {
-                storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::MarkovChainSettings>().getEquationSolver();
-                switch (equationSolver) {
-                    case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::EliminationLinearEquationSolver<storm::RationalNumber>>(matrix);
-                    default: return std::make_unique<storm::solver::EigenLinearEquationSolver<storm::RationalNumber>>(matrix);
-                }
-            }
-            
-            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> GeneralLinearEquationSolverFactory<storm::RationalFunction>::create(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix) const {
-                return selectSolver(matrix);
-            }
-
-            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> GeneralLinearEquationSolverFactory<storm::RationalFunction>::create(storm::storage::SparseMatrix<storm::RationalFunction>&& matrix) const {
-                return selectSolver(std::move(matrix));
-            }
-
-            template<typename MatrixType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> GeneralLinearEquationSolverFactory<storm::RationalFunction>::selectSolver(MatrixType&& matrix) const {
-                storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::MarkovChainSettings>().getEquationSolver();
-                switch (equationSolver) {
-                    case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::EliminationLinearEquationSolver<storm::RationalFunction>>(matrix);
-                    default: return std::make_unique<storm::solver::EigenLinearEquationSolver<storm::RationalFunction>>(matrix);
-                }
-            }
-            
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
-                return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(matrix);
-            }
-
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
-                return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(std::move(matrix));
-            }
-
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EigenLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
-                return std::make_unique<storm::solver::EigenLinearEquationSolver<ValueType>>(matrix);
-            }
-
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EigenLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
-                return std::make_unique<storm::solver::EigenLinearEquationSolver<ValueType>>(std::move(matrix));
-            }
-            
-            template<typename ValueType>
-            NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory() {
-                switch (storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getLinearEquationSystemMethod()) {
-                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Jacobi:
-                        this->method = storm::solver::NativeLinearEquationSolverSolutionMethod::Jacobi;
-                        break;
-                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel:
-                        this->method = storm::solver::NativeLinearEquationSolverSolutionMethod::GaussSeidel;
-                        break;
-                    case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::SOR:
-                        this->method = storm::solver::NativeLinearEquationSolverSolutionMethod::SOR;
-                        break;
-                }
-            }
-            
-            template<typename ValueType>
-            NativeLinearEquationSolverFactory<ValueType>::NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolverSolutionMethod method) : method(method) {
-                // Intentionally left empty.
-            }
-            
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
-                return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(matrix, method);
-            }
-
-            template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType>&& matrix) const {
-                return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(std::move(matrix), method);
-            }
-            
             template<typename ValueType>
             MinMaxLinearEquationSolverFactory<ValueType>::MinMaxLinearEquationSolverFactory(storm::solver::EquationSolverTypeSelection solver)
             {
@@ -269,21 +157,8 @@ namespace storm {
             template class SymbolicMinMaxLinearEquationSolverFactory<storm::dd::DdType::Sylvan, double>;
             template class SymbolicGameSolverFactory<storm::dd::DdType::CUDD, double>;
             template class SymbolicGameSolverFactory<storm::dd::DdType::Sylvan, double>;
-            template class LinearEquationSolverFactory<double>;
-            template class GeneralLinearEquationSolverFactory<double>;
-            template class GmmxxLinearEquationSolverFactory<double>;
-            template class EigenLinearEquationSolverFactory<double>;
-            template class NativeLinearEquationSolverFactory<double>;
             template class MinMaxLinearEquationSolverFactory<double>;
             template class GameSolverFactory<double>;
-
-            template class LinearEquationSolverFactory<storm::RationalNumber>;
-            template class GeneralLinearEquationSolverFactory<storm::RationalNumber>;
-            template class EigenLinearEquationSolverFactory<storm::RationalNumber>;
-
-            template class LinearEquationSolverFactory<storm::RationalFunction>;
-            template class GeneralLinearEquationSolverFactory<storm::RationalFunction>;
-            template class EigenLinearEquationSolverFactory<storm::RationalFunction>;
         }
     }
 }
diff --git a/src/utility/solver.h b/src/utility/solver.h
index a5fab02c7..d47c304a7 100644
--- a/src/utility/solver.h
+++ b/src/utility/solver.h
@@ -33,11 +33,6 @@ namespace storm {
         
         class LpSolver;
         class SmtSolver;
-        
-        template<typename ValueType>
-        class NativeLinearEquationSolver;
-        
-        enum class NativeLinearEquationSolverSolutionMethod;
     }
 
     namespace storage {
@@ -77,87 +72,6 @@ namespace storm {
             public:
                 virtual std::unique_ptr<storm::solver::SymbolicGameSolver<Type, ValueType>> create(storm::dd::Add<Type, ValueType> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::set<storm::expressions::Variable> const& player1Variables, std::set<storm::expressions::Variable> const& player2Variables) const;
             };
-
-            template<typename ValueType>
-            class LinearEquationSolverFactory {
-            public:
-                /*!
-                 * Creates a new linear equation solver instance with the given matrix.
-                 *
-                 * @param matrix The matrix that defines the equation system.
-                 * @return A pointer to the newly created solver.
-                 */
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const = 0;
-                
-                /*!
-                 * Creates a new linear equation solver instance with the given matrix. The caller gives up posession of the
-                 * matrix by calling this function.
-                 *
-                 * @param matrix The matrix that defines the equation system.
-                 * @return A pointer to the newly created solver.
-                 */
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const;
-            };
-            
-            template<typename ValueType>
-            class GeneralLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
-            public:
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
-                
-            private:
-                template<typename MatrixType>
-                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> selectSolver(MatrixType&& matrix) const;
-            };
-
-            template<>
-            class GeneralLinearEquationSolverFactory<storm::RationalNumber> : public LinearEquationSolverFactory<storm::RationalNumber> {
-            public:
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> create(storm::storage::SparseMatrix<storm::RationalNumber> const& matrix) const override;
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> create(storm::storage::SparseMatrix<storm::RationalNumber>&& matrix) const override;
-                
-            private:
-                template<typename MatrixType>
-                std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalNumber>> selectSolver(MatrixType&& matrix) const;
-            };
-
-            template<>
-            class GeneralLinearEquationSolverFactory<storm::RationalFunction> : public LinearEquationSolverFactory<storm::RationalFunction> {
-            public:
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> create(storm::storage::SparseMatrix<storm::RationalFunction> const& matrix) const override;
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> create(storm::storage::SparseMatrix<storm::RationalFunction>&& matrix) const override;
-
-            private:
-                template<typename MatrixType>
-                std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> selectSolver(MatrixType&& matrix) const;
-            };
-            
-            template<typename ValueType>
-            class NativeLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
-            public:
-                NativeLinearEquationSolverFactory();
-                NativeLinearEquationSolverFactory(storm::solver::NativeLinearEquationSolverSolutionMethod method);
-                
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
-                
-            private:
-                storm::solver::NativeLinearEquationSolverSolutionMethod method;
-            };
-            
-            template<typename ValueType>
-            class GmmxxLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
-            public:
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
-            };
-            
-            template<typename ValueType>
-            class EigenLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
-            public:
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
-                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType>&& matrix) const override;
-            };
             
             template<typename ValueType>
             class MinMaxLinearEquationSolverFactory {
diff --git a/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp
index 113d0b1ba..8caa788bb 100644
--- a/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/EigenDtmcPrctlModelCheckerTest.cpp
@@ -3,7 +3,7 @@
 
 #include "src/parser/FormulaParser.h"
 #include "src/logic/Formulas.h"
-#include "src/utility/solver.h"
+#include "src/solver/EigenLinearEquationSolver.h"
 #include "src/models/sparse/StandardRewardModel.h"
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
@@ -29,7 +29,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, Die) {
     ASSERT_EQ(dtmc->getNumberOfStates(), 13ull);
     ASSERT_EQ(dtmc->getNumberOfTransitions(), 20ull);
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
     
@@ -77,7 +77,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, Die_RationalNumber) {
     ASSERT_EQ(dtmc->getNumberOfStates(), 13ull);
     ASSERT_EQ(dtmc->getNumberOfTransitions(), 20ull);
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<storm::RationalNumber>> checker(*dtmc, std::make_unique<storm::utility::solver::EigenLinearEquationSolverFactory<storm::RationalNumber>>());
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<storm::RationalNumber>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<storm::RationalNumber>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
     
@@ -129,7 +129,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, Die_RationalFunction) {
     ASSERT_EQ(variables.size(), 1ull);
     instantiation.emplace(*variables.begin(), storm::RationalNumber(1) / storm::RationalNumber(2));
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>> checker(*dtmc, std::make_unique<storm::utility::solver::EigenLinearEquationSolverFactory<storm::RationalFunction>>());
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<storm::RationalFunction>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<storm::RationalFunction>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
     
@@ -174,7 +174,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, Crowds) {
     ASSERT_EQ(8607ull, dtmc->getNumberOfStates());
     ASSERT_EQ(15113ull, dtmc->getNumberOfTransitions());
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::solver::LinearEquationSolverFactory<double>>(new storm::solver::EigenLinearEquationSolverFactory<double>()));
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]");
     
@@ -211,7 +211,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, SynchronousLeader) {
     ASSERT_EQ(12400ull, dtmc->getNumberOfStates());
     ASSERT_EQ(16495ull, dtmc->getNumberOfTransitions());
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]");
     
@@ -254,7 +254,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, LRASingleBscc) {
         
         dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
         
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<double>>());
         
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
         
@@ -278,7 +278,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, LRASingleBscc) {
         
         dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
         
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<double>>());
         
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
         
@@ -302,7 +302,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, LRASingleBscc) {
         
         dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
         
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<double>>());
         
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
         
@@ -365,7 +365,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, LRA) {
         
         dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
         
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<double>>());
         
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
         
@@ -394,7 +394,7 @@ TEST(EigenDtmcPrctlModelCheckerTest, Conditional) {
     
     std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc = model->as<storm::models::sparse::Dtmc<double>>();
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::EigenLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::EigenLinearEquationSolverFactory<double>>());
     
     // A parser that we use for conveniently constructing the formulas.
     storm::parser::FormulaParser formulaParser;
diff --git a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
index 25763b37f..8f6814ded 100644
--- a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
@@ -35,7 +35,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Cluster) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=100 !\"minimum\"]");
@@ -111,7 +111,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Embedded) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10000 \"down\"]");
@@ -173,7 +173,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Polling) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=?[ F<=10 \"target\"]");
@@ -214,7 +214,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Tandem) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10 \"network_full\" ]");
diff --git a/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
index bb96f4511..2290a1cec 100644
--- a/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
@@ -3,7 +3,7 @@
 
 #include "src/parser/FormulaParser.h"
 #include "src/logic/Formulas.h"
-#include "src/utility/solver.h"
+#include "src/solver/GmmxxLinearEquationSolver.h"
 #include "src/models/sparse/StandardRewardModel.h"
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
@@ -29,7 +29,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, Die) {
     ASSERT_EQ(dtmc->getNumberOfStates(), 13ull);
     ASSERT_EQ(dtmc->getNumberOfTransitions(), 20ull);
 
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
 
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
 
@@ -73,7 +73,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, Crowds) {
     ASSERT_EQ(8607ull, dtmc->getNumberOfStates());
     ASSERT_EQ(15113ull, dtmc->getNumberOfTransitions());
 
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
 
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]");
 
@@ -110,7 +110,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, SynchronousLeader) {
     ASSERT_EQ(12400ull, dtmc->getNumberOfStates());
     ASSERT_EQ(16495ull, dtmc->getNumberOfTransitions());
 
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
 
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]");
 
@@ -153,7 +153,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
         dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
 
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
 
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
 
@@ -177,7 +177,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
         dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
 
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
 
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
 
@@ -201,7 +201,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
         dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
 
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
 
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
 
@@ -264,7 +264,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, LRA) {
 
         dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
 
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
 
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
 
@@ -291,7 +291,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, Conditional) {
 
     std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc = model->as<storm::models::sparse::Dtmc<double>>();
 
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
 
     // A parser that we use for conveniently constructing the formulas.
     storm::parser::FormulaParser formulaParser;
diff --git a/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp
index 52196371d..b5e26c9ca 100644
--- a/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxHybridCtmcCslModelCheckerTest.cpp
@@ -7,7 +7,7 @@
 #include "src/builder/DdPrismModelBuilder.h"
 #include "src/storage/dd/DdType.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/GmmxxLinearEquationSolver.h"
 #include "src/models/symbolic/StandardRewardModel.h"
 #include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
 #include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
@@ -43,7 +43,7 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Cluster_Cudd) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=100 !\"minimum\"]");
@@ -140,7 +140,7 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Cluster_Sylvan) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=100 !\"minimum\"]");
@@ -237,7 +237,7 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Embedded_Cudd) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10000 \"down\"]");
@@ -316,7 +316,7 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Embedded_Sylvan) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10000 \"down\"]");
@@ -388,7 +388,7 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Polling_Cudd) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=?[ F<=10 \"target\"]");
@@ -424,7 +424,7 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Polling_Sylvan) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=?[ F<=10 \"target\"]");
@@ -474,7 +474,7 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Tandem_Cudd) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10 \"network_full\" ]");
@@ -562,7 +562,7 @@ TEST(GmmxxHybridCtmcCslModelCheckerTest, Tandem_Sylvan) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10 \"network_full\" ]");
diff --git a/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp
index 0c4679382..a8a0857ad 100644
--- a/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxHybridDtmcPrctlModelCheckerTest.cpp
@@ -3,7 +3,7 @@
 
 #include "src/parser/FormulaParser.h"
 #include "src/logic/Formulas.h"
-#include "src/utility/solver.h"
+#include "src/solver/GmmxxLinearEquationSolver.h"
 #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
 #include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
 #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
@@ -40,7 +40,7 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, Die_Cudd) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
     
@@ -101,7 +101,7 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, Die_Sylvan) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
     
@@ -154,7 +154,7 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, Crowds_Cudd) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]");
     
@@ -198,7 +198,7 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, Crowds_Sylvan) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]");
     
@@ -250,7 +250,7 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, SynchronousLeader_Cudd) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]");
     
@@ -302,7 +302,7 @@ TEST(GmmxxHybridDtmcPrctlModelCheckerTest, SynchronousLeader_Sylvan) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::make_unique<storm::solver::GmmxxLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]");
     
diff --git a/test/functional/modelchecker/NativeCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/NativeCtmcCslModelCheckerTest.cpp
index 65d58c19f..29a937a5e 100644
--- a/test/functional/modelchecker/NativeCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/NativeCtmcCslModelCheckerTest.cpp
@@ -6,7 +6,7 @@
 #include "src/logic/Formulas.h"
 #include "src/builder/ExplicitModelBuilder.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/NativeLinearEquationSolver.h"
 #include "src/models/sparse/StandardRewardModel.h"
 #include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
@@ -32,7 +32,7 @@ TEST(NativeCtmcCslModelCheckerTest, Cluster) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
 
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=100 !\"minimum\"]");
@@ -103,7 +103,7 @@ TEST(NativeCtmcCslModelCheckerTest, Embedded) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
 
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10000 \"down\"]");
@@ -158,7 +158,7 @@ TEST(NativeCtmcCslModelCheckerTest, Polling) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=?[ F<=10 \"target\"]");
@@ -192,7 +192,7 @@ TEST(NativeCtmcCslModelCheckerTest, Tandem) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<storm::models::sparse::Ctmc<double>> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10 \"network_full\" ]");
diff --git a/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
index b53c080e2..b8b5bda5a 100644
--- a/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
@@ -4,7 +4,7 @@
 #include "src/parser/FormulaParser.h"
 #include "src/settings/SettingMemento.h"
 #include "src/logic/Formulas.h"
-#include "src/utility/solver.h"
+#include "src/solver/NativeLinearEquationSolver.h"
 #include "src/models/sparse/StandardRewardModel.h"
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
@@ -29,7 +29,7 @@ TEST(NativeDtmcPrctlModelCheckerTest, Die) {
     ASSERT_EQ(dtmc->getNumberOfStates(), 13ull);
     ASSERT_EQ(dtmc->getNumberOfTransitions(), 20ull);
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
     
@@ -73,7 +73,7 @@ TEST(NativeDtmcPrctlModelCheckerTest, Crowds) {
     ASSERT_EQ(8607ull, dtmc->getNumberOfStates());
     ASSERT_EQ(15113ull, dtmc->getNumberOfTransitions());
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]");
     
@@ -110,7 +110,7 @@ TEST(NativeDtmcPrctlModelCheckerTest, SynchronousLeader) {
     ASSERT_EQ(12400ull, dtmc->getNumberOfStates());
     ASSERT_EQ(16495ull, dtmc->getNumberOfTransitions());
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]");
     
@@ -153,7 +153,10 @@ TEST(NativeDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
 
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolverSolutionMethod::SOR, 0.9)));
+        auto factory = std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>();
+        factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings<double>::SolutionMethod::SOR);
+        factory->getSettings().setOmega(0.9);
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::move(factory));
 
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
         
@@ -177,7 +180,10 @@ TEST(NativeDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
 
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolverSolutionMethod::SOR, 0.9)));
+        auto factory = std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>();
+        factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings<double>::SolutionMethod::SOR);
+        factory->getSettings().setOmega(0.9);
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::move(factory));
 
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
         
@@ -201,7 +207,10 @@ TEST(NativeDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
 		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
 
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolverSolutionMethod::SOR, 0.9)));
+        auto factory = std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>();
+        factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings<double>::SolutionMethod::SOR);
+        factory->getSettings().setOmega(0.9);
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::move(factory));
 
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
         
@@ -216,7 +225,7 @@ TEST(NativeDtmcPrctlModelCheckerTest, LRASingleBscc) {
 
 TEST(NativeDtmcPrctlModelCheckerTest, LRA) {
 	storm::storage::SparseMatrixBuilder<double> matrixBuilder;
-	std::shared_ptr<storm::models::sparse::Dtmc<double>> mdp;
+    std::shared_ptr<storm::models::sparse::Dtmc<double>> dtmc;
 
     // A parser that we use for conveniently constructing the formulas.
     storm::parser::FormulaParser formulaParser;
@@ -262,9 +271,12 @@ TEST(NativeDtmcPrctlModelCheckerTest, LRA) {
 		ap.addLabelToState("a", 13);
 		ap.addLabelToState("a", 14);
 
-		mdp.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
+		dtmc.reset(new storm::models::sparse::Dtmc<double>(transitionMatrix, ap));
 
-        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*mdp, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>(storm::solver::NativeLinearEquationSolverSolutionMethod::SOR, 0.9)));
+        auto factory = std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>();
+        factory->getSettings().setSolutionMethod(storm::solver::NativeLinearEquationSolverSettings<double>::SolutionMethod::SOR);
+        factory->getSettings().setOmega(0.9);
+        storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::move(factory));
 
         std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("LRA=? [\"a\"]");
         
diff --git a/test/functional/modelchecker/NativeHybridCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/NativeHybridCtmcCslModelCheckerTest.cpp
index d2e3ff7d6..58977dc60 100644
--- a/test/functional/modelchecker/NativeHybridCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/NativeHybridCtmcCslModelCheckerTest.cpp
@@ -7,7 +7,7 @@
 #include "src/builder/DdPrismModelBuilder.h"
 #include "src/storage/dd/DdType.h"
 
-#include "src/utility/solver.h"
+#include "src/solver/NativeLinearEquationSolver.h"
 #include "src/models/symbolic/StandardRewardModel.h"
 #include "src/modelchecker/csl/HybridCtmcCslModelChecker.h"
 #include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
@@ -42,7 +42,7 @@ TEST(NativeHybridCtmcCslModelCheckerTest, Cluster_Cudd) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=100 !\"minimum\"]");
@@ -139,7 +139,7 @@ TEST(NativeHybridCtmcCslModelCheckerTest, Cluster_Sylvan) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=100 !\"minimum\"]");
@@ -236,7 +236,7 @@ TEST(NativeHybridCtmcCslModelCheckerTest, Embedded_Cudd) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10000 \"down\"]");
@@ -315,7 +315,7 @@ TEST(NativeHybridCtmcCslModelCheckerTest, Embedded_Sylvan) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10000 \"down\"]");
@@ -387,7 +387,7 @@ TEST(NativeHybridCtmcCslModelCheckerTest, Polling_Cudd) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=?[ F<=10 \"target\"]");
@@ -423,7 +423,7 @@ TEST(NativeHybridCtmcCslModelCheckerTest, Polling_Sylvan) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=?[ F<=10 \"target\"]");
@@ -473,7 +473,7 @@ TEST(NativeHybridCtmcCslModelCheckerTest, Tandem_Cudd) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::CUDD>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::CUDD, double> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10 \"network_full\" ]");
@@ -563,7 +563,7 @@ TEST(NativeHybridCtmcCslModelCheckerTest, Tandem_Sylvan) {
     std::shared_ptr<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>> ctmc = model->as<storm::models::symbolic::Ctmc<storm::dd::DdType::Sylvan>>();
     
     // Create model checker.
-    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridCtmcCslModelChecker<storm::dd::DdType::Sylvan, double> modelchecker(*ctmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     // Start checking properties.
     formula = formulaParser.parseSingleFormulaFromString("P=? [ F<=10 \"network_full\" ]");
diff --git a/test/functional/modelchecker/NativeHybridDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeHybridDtmcPrctlModelCheckerTest.cpp
index 072745694..bac3f8aa3 100644
--- a/test/functional/modelchecker/NativeHybridDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/NativeHybridDtmcPrctlModelCheckerTest.cpp
@@ -3,7 +3,7 @@
 
 #include "src/parser/FormulaParser.h"
 #include "src/logic/Formulas.h"
-#include "src/utility/solver.h"
+#include "src/solver/NativeLinearEquationSolver.h"
 #include "src/modelchecker/prctl/HybridDtmcPrctlModelChecker.h"
 #include "src/modelchecker/results/HybridQuantitativeCheckResult.h"
 #include "src/modelchecker/results/SymbolicQualitativeCheckResult.h"
@@ -41,7 +41,7 @@ TEST(NativeHybridDtmcPrctlModelCheckerTest, Die_Cudd) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
     
@@ -102,7 +102,7 @@ TEST(NativeHybridDtmcPrctlModelCheckerTest, Die_Sylvan) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"one\"]");
     
@@ -155,7 +155,7 @@ TEST(NativeHybridDtmcPrctlModelCheckerTest, Crowds_Cudd) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]");
     
@@ -199,7 +199,7 @@ TEST(NativeHybridDtmcPrctlModelCheckerTest, Crowds_Sylvan) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"observe0Greater1\"]");
     
@@ -251,7 +251,7 @@ TEST(NativeHybridDtmcPrctlModelCheckerTest, SynchronousLeader_Cudd) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::CUDD>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::CUDD, double> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]");
     
@@ -303,7 +303,7 @@ TEST(NativeHybridDtmcPrctlModelCheckerTest, SynchronousLeader_Sylvan) {
     
     std::shared_ptr<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>> dtmc = model->as<storm::models::symbolic::Dtmc<storm::dd::DdType::Sylvan>>();
     
-    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::HybridDtmcPrctlModelChecker<storm::dd::DdType::Sylvan, double> checker(*dtmc, std::make_unique<storm::solver::NativeLinearEquationSolverFactory<double>>());
     
     std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("P=? [F \"elected\"]");
     
diff --git a/test/functional/solver/EigenLinearEquationSolverTest.cpp b/test/functional/solver/EigenLinearEquationSolverTest.cpp
index 65710ca1b..2bf575c2a 100644
--- a/test/functional/solver/EigenLinearEquationSolverTest.cpp
+++ b/test/functional/solver/EigenLinearEquationSolverTest.cpp
@@ -54,9 +54,10 @@ TEST(EigenLinearEquationSolver, SparseLU) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::SparseLU, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None));
-    
-    storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::SparseLU, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None);
+    storm::solver::EigenLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::EigenLinearEquationSolverSettings<double>::SolutionMethod::SparseLU);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::EigenLinearEquationSolverSettings<double>::Preconditioner::None);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
@@ -82,9 +83,7 @@ TEST(EigenLinearEquationSolver, SparseLU_RationalNumber) {
     std::vector<storm::RationalNumber> x(3);
     std::vector<storm::RationalNumber> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<storm::RationalNumber> solver(A, storm::solver::EigenLinearEquationSolver<storm::RationalNumber>::SolutionMethod::SparseLU));
-    
-    storm::solver::EigenLinearEquationSolver<storm::RationalNumber> solver(A, storm::solver::EigenLinearEquationSolver<storm::RationalNumber>::SolutionMethod::SparseLU);
+    storm::solver::EigenLinearEquationSolver<storm::RationalNumber> solver(A);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_TRUE(storm::utility::isOne(x[0]));
     ASSERT_TRUE(x[1] == 3);
@@ -110,9 +109,7 @@ TEST(EigenLinearEquationSolver, SparseLU_RationalFunction) {
     std::vector<storm::RationalFunction> x(3);
     std::vector<storm::RationalFunction> b = {storm::RationalFunction(16), storm::RationalFunction(-4), storm::RationalFunction(-7)};
     
-    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<storm::RationalFunction> solver(A, storm::solver::EigenLinearEquationSolver<storm::RationalFunction>::SolutionMethod::SparseLU));
-    
-    storm::solver::EigenLinearEquationSolver<storm::RationalFunction> solver(A, storm::solver::EigenLinearEquationSolver<storm::RationalFunction>::SolutionMethod::SparseLU);
+    storm::solver::EigenLinearEquationSolver<storm::RationalFunction> solver(A);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_TRUE(storm::utility::isOne(x[0]));
     ASSERT_TRUE(x[1] == storm::RationalFunction(3));
@@ -138,9 +135,12 @@ TEST(DISABLED_EigenLinearEquationSolver, BiCGSTAB) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None));
-    
-    storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::None);
+    storm::solver::EigenLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::EigenLinearEquationSolverSettings<double>::SolutionMethod::Bicgstab);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::EigenLinearEquationSolverSettings<double>::Preconditioner::None);
+
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
@@ -166,9 +166,11 @@ TEST(DISABLED_EigenLinearEquationSolver, BiCGSTAB_Ilu) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Ilu));
-    
-    storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Ilu);
+    storm::solver::EigenLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::EigenLinearEquationSolverSettings<double>::SolutionMethod::Bicgstab);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::EigenLinearEquationSolverSettings<double>::Preconditioner::Ilu);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
@@ -194,9 +196,11 @@ TEST(DISABLED_EigenLinearEquationSolver, BiCGSTAB_Diagonal) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Diagonal));
-    
-    storm::solver::EigenLinearEquationSolver<double> solver(A, storm::solver::EigenLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::EigenLinearEquationSolver<double>::Preconditioner::Diagonal);
+    storm::solver::EigenLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::EigenLinearEquationSolverSettings<double>::SolutionMethod::Bicgstab);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::EigenLinearEquationSolverSettings<double>::Preconditioner::Diagonal);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::EigenEquationSolverSettings>().getPrecision());
diff --git a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
index 2820678d2..181d3429a 100644
--- a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
+++ b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
@@ -53,9 +53,13 @@ TEST(GmmxxLinearEquationSolver, gmres) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50));
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings<double>::SolutionMethod::Gmres);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<double>::Preconditioner::None);
+    solver.getSettings().setNumberOfIterationsUntilRestart(50);
     
-    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
@@ -81,10 +85,19 @@ TEST(GmmxxLinearEquationSolver, qmr) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Qmr, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None));
-    
-    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Qmr, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50);
-    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings<double>::SolutionMethod::Qmr);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<double>::Preconditioner::None);
+    
+    storm::solver::GmmxxLinearEquationSolver<double> solver2(A);
+    solver2.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings<double>::SolutionMethod::Qmr);
+    solver2.getSettings().setPrecision(1e-6);
+    solver2.getSettings().setMaximalNumberOfIterations(10000);
+    solver2.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<double>::Preconditioner::None);
+    
+    ASSERT_NO_THROW(solver2.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
@@ -109,9 +122,12 @@ TEST(GmmxxLinearEquationSolver, bicgstab) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None));
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings<double>::SolutionMethod::Bicgstab);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<double>::Preconditioner::None);
     
-    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
@@ -137,9 +153,11 @@ TEST(GmmxxLinearEquationSolver, jacobi) {
     std::vector<double> x(3);
     std::vector<double> b = {11, -16, 1};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Jacobi, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None));
-    
-    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Jacobi, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None);
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings<double>::SolutionMethod::Jacobi);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<double>::Preconditioner::None);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
@@ -165,9 +183,12 @@ TEST(GmmxxLinearEquationSolver, gmresilu) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Ilu, true, 50));
-    
-    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Ilu, true, 50);
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings<double>::SolutionMethod::Gmres);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<double>::Preconditioner::Ilu);
+    solver.getSettings().setNumberOfIterationsUntilRestart(50);
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
@@ -193,9 +214,13 @@ TEST(GmmxxLinearEquationSolver, gmresdiag) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Diagonal, true, 50));
-    
-    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Diagonal, true, 50);
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A);
+    solver.getSettings().setSolutionMethod(storm::solver::GmmxxLinearEquationSolverSettings<double>::SolutionMethod::Gmres);
+    solver.getSettings().setPrecision(1e-6);
+    solver.getSettings().setMaximalNumberOfIterations(10000);
+    solver.getSettings().setPreconditioner(storm::solver::GmmxxLinearEquationSolverSettings<double>::Preconditioner::Diagonal);
+    solver.getSettings().setNumberOfIterationsUntilRestart(50);
+
     ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
diff --git a/test/performance/modelchecker/GmmxxDtmcPrctModelCheckerTest.cpp b/test/performance/modelchecker/GmmxxDtmcPrctModelCheckerTest.cpp
index 8f3f9607f..18d0f24bd 100644
--- a/test/performance/modelchecker/GmmxxDtmcPrctModelCheckerTest.cpp
+++ b/test/performance/modelchecker/GmmxxDtmcPrctModelCheckerTest.cpp
@@ -21,7 +21,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, Crowds) {
 	ASSERT_EQ(2036647ull, dtmc->getNumberOfStates());
 	ASSERT_EQ(7362293ull, dtmc->getNumberOfTransitions());
 
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
 
     // A parser that we use for conveniently constructing the formulas.
     storm::parser::FormulaParser formulaParser;
@@ -59,7 +59,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, SynchronousLeader) {
 	ASSERT_EQ(1312334ull, dtmc->getNumberOfStates());
 	ASSERT_EQ(1574477ull, dtmc->getNumberOfTransitions());
 
-	storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
+	storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
 
     // A parser that we use for conveniently constructing the formulas.
     storm::parser::FormulaParser formulaParser;
diff --git a/test/performance/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp b/test/performance/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
index 94635ee95..95a112a56 100644
--- a/test/performance/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
+++ b/test/performance/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp
@@ -21,7 +21,7 @@ TEST(NativeDtmcPrctlModelCheckerTest, Crowds) {
     ASSERT_EQ(2036647ull, dtmc->getNumberOfStates());
     ASSERT_EQ(7362293ull, dtmc->getNumberOfTransitions());
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
     
     // A parser that we use for conveniently constructing the formulas.
     storm::parser::FormulaParser formulaParser;
@@ -59,7 +59,7 @@ TEST(NativeDtmcPrctlModelCheckerTest, SynchronousLeader) {
     ASSERT_EQ(1312334ull, dtmc->getNumberOfStates());
     ASSERT_EQ(1574477ull, dtmc->getNumberOfTransitions());
     
-    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<storm::models::sparse::Dtmc<double>> checker(*dtmc, std::unique_ptr<storm::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
 
     // A parser that we use for conveniently constructing the formulas.
     storm::parser::FormulaParser formulaParser;