diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
index ed78b8669..80d1147e6 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
@@ -19,18 +19,17 @@
 namespace storm {
     namespace modelchecker {
         template<class ValueType>
-        SparseCtmcCslModelChecker<ValueType>::SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolver(storm::utility::solver::getLinearEquationSolver<ValueType>()) {
+        SparseCtmcCslModelChecker<ValueType>::SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolverFactory(new storm::utility::solver::LinearEquationSolverFactory<ValueType>()) {
             // Intentionally left empty.
         }
         
         template<class ValueType>
-        SparseCtmcCslModelChecker<ValueType>::SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>&& linearEquationSolver) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolver(std::move(linearEquationSolver)) {
+        SparseCtmcCslModelChecker<ValueType>::SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
         template<class ValueType>
         bool SparseCtmcCslModelChecker<ValueType>::canHandle(storm::logic::Formula const& formula) const {
-            // FIXME: refine.
             return formula.isCslStateFormula() || formula.isCslPathFormula() || formula.isRewardPathFormula();
         }
         
@@ -59,7 +58,7 @@ namespace storm {
         std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            std::vector<ValueType> result = SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilitiesHelper(this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), subResult.getTruthValuesVector(), *this->linearEquationSolver);
+            std::vector<ValueType> result = SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilitiesHelper(this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory);
             return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(std::move(result)));
         }
         
@@ -69,7 +68,7 @@ namespace storm {
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
             ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolver)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory)));
         }
         
         template<class ValueType>
@@ -82,7 +81,7 @@ namespace storm {
             // If the time bounds are [0, inf], we rather call untimed reachability.
             storm::utility::ConstantsComparator<ValueType> comparator;
             if (comparator.isZero(lowerBound) && comparator.isInfinity(upperBound)) {
-                return this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, qualitative, *this->linearEquationSolver);
+                return this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, qualitative, *this->linearEquationSolverFactory);
             }
             
             // From this point on, we know that we have to solve a more complicated problem [t, t'] with either t != 0
@@ -124,7 +123,7 @@ namespace storm {
                         std::vector<ValueType> psiStateValues(statesWithProbabilityGreater0.getNumberOfSetBits(), storm::utility::zero<ValueType>());
                         storm::utility::vector::setVectorValues(psiStateValues, psiStates % statesWithProbabilityGreater0, storm::utility::one<ValueType>());
                         
-                        std::vector<ValueType> subresult = this->computeTransientProbabilities(uniformizedMatrix, upperBound, uniformizationRate, psiStateValues, *this->linearEquationSolver);
+                        std::vector<ValueType> subresult = this->computeTransientProbabilities(uniformizedMatrix, upperBound, uniformizationRate, psiStateValues, *this->linearEquationSolverFactory);
                         result = std::vector<ValueType>(this->getModel().getNumberOfStates(), storm::utility::zero<ValueType>());
                         storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, subresult);
                     } else if (comparator.isInfinity(upperBound)) {
@@ -132,7 +131,7 @@ namespace storm {
                         
                         // Start by computing the (unbounded) reachability probabilities of reaching psi states while
                         // staying in phi states.
-                        result = this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), backwardTransitions, phiStates, psiStates, qualitative, *this->linearEquationSolver);
+                        result = this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), backwardTransitions, phiStates, psiStates, qualitative, *this->linearEquationSolverFactory);
                         
                         // Determine the set of states that must be considered further.
                         storm::storage::BitVector relevantStates = storm::utility::vector::filterGreaterZero(result);
@@ -151,7 +150,7 @@ namespace storm {
                         storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), relevantStates, storm::storage::BitVector(result.size()), uniformizationRate, exitRates);
                         
                         // Compute the transient probabilities.
-                        subResult = this->computeTransientProbabilities(uniformizedMatrix, lowerBound, uniformizationRate, subResult, *this->linearEquationSolver);
+                        subResult = this->computeTransientProbabilities(uniformizedMatrix, lowerBound, uniformizationRate, subResult, *this->linearEquationSolverFactory);
                         
                         // Fill in the correct values.
                         storm::utility::vector::setVectorValues(result, ~relevantStates, storm::utility::zero<ValueType>());
@@ -173,7 +172,7 @@ namespace storm {
                         // Start by computing the transient probabilities of reaching a psi state in time t' - t.
                         std::vector<ValueType> psiStateValues(statesWithProbabilityGreater0.getNumberOfSetBits(), storm::utility::zero<ValueType>());
                         storm::utility::vector::setVectorValues(psiStateValues, psiStates % statesWithProbabilityGreater0, storm::utility::one<ValueType>());
-                        std::vector<ValueType> subresult = this->computeTransientProbabilities(uniformizedMatrix, upperBound - lowerBound, uniformizationRate, psiStateValues, *this->linearEquationSolver);
+                        std::vector<ValueType> subresult = this->computeTransientProbabilities(uniformizedMatrix, upperBound - lowerBound, uniformizationRate, psiStateValues, *this->linearEquationSolverFactory);
                         
                         // Determine the set of states that must be considered further.
                         storm::storage::BitVector relevantStates = storm::utility::vector::filterGreaterZero(subresult);
@@ -193,7 +192,7 @@ namespace storm {
                         
                         // Finally, we compute the second set of transient probabilities.
                         uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), relevantStates, storm::storage::BitVector(this->getModel().getNumberOfStates()), uniformizationRate, exitRates);
-                        newSubresult = this->computeTransientProbabilities(uniformizedMatrix, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolver);
+                        newSubresult = this->computeTransientProbabilities(uniformizedMatrix, lowerBound, uniformizationRate, newSubresult, *this->linearEquationSolverFactory);
                         
                         // Fill in the correct values.
                         result = std::vector<ValueType>(this->getModel().getNumberOfStates(), storm::utility::zero<ValueType>());
@@ -241,7 +240,7 @@ namespace storm {
         
         template<class ValueType>
         template<bool computeCumulativeReward>
-        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver) const {
+        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) const {
             ValueType lambda = timeBound * uniformizationRate;
             
             // If no time can pass, the current values are the result.
@@ -285,15 +284,17 @@ namespace storm {
             }
             std::vector<ValueType> multiplicationResult(result.size());
             
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(uniformizedMatrix);
+            
             if (!computeCumulativeReward && std::get<0>(foxGlynnResult) > 1) {
                 // Perform the matrix-vector multiplications (without adding).
-                linearEquationSolver.performMatrixVectorMultiplication(uniformizedMatrix, values, nullptr, std::get<0>(foxGlynnResult) - 1, &multiplicationResult);
+                solver->performMatrixVectorMultiplication(values, nullptr, std::get<0>(foxGlynnResult) - 1, &multiplicationResult);
             } else if (computeCumulativeReward) {
                 std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; };
                 
                 // For the iterations below the left truncation point, we need to add and scale the result with the uniformization rate.
                 for (uint_fast64_t index = 1; index < startingIteration; ++index) {
-                    linearEquationSolver.performMatrixVectorMultiplication(uniformizedMatrix, values, nullptr, 1, &multiplicationResult);
+                    solver->performMatrixVectorMultiplication(values, nullptr, 1, &multiplicationResult);
                     storm::utility::vector::applyPointwiseInPlace(result, values, addAndScale);
                 }
             }
@@ -303,7 +304,7 @@ namespace storm {
             ValueType weight = 0;
             std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; };
             for (uint_fast64_t index = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) {
-                linearEquationSolver.performMatrixVectorMultiplication(uniformizedMatrix, values, nullptr, 1, &multiplicationResult);
+                solver->performMatrixVectorMultiplication(values, nullptr, 1, &multiplicationResult);
                 
                 weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
                 storm::utility::vector::applyPointwiseInPlace(result, values, addAndScale);
@@ -313,8 +314,8 @@ namespace storm {
         }
         
         template<class ValueType>
-        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeUntilProbabilitiesHelper(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::LinearEquationSolver<ValueType> const& linearEquationSolver) {
-            return SparseDtmcPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(transitionMatrix, backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolver);
+        std::vector<ValueType> SparseCtmcCslModelChecker<ValueType>::computeUntilProbabilitiesHelper(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) {
+            return SparseDtmcPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(transitionMatrix, backwardTransitions, phiStates, psiStates, qualitative, linearEquationSolverFactory);
         }
         
         template<class ValueType>
@@ -340,7 +341,7 @@ namespace storm {
                 STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
                 
                 storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), storm::storage::BitVector(this->getModel().getNumberOfStates()), uniformizationRate, this->getModel().getExitRateVector());
-                result = this->computeTransientProbabilities(uniformizedMatrix, timeBound, uniformizationRate, result, *this->linearEquationSolver);
+                result = this->computeTransientProbabilities(uniformizedMatrix, timeBound, uniformizationRate, result, *this->linearEquationSolverFactory);
             }
             
             return result;
@@ -385,7 +386,7 @@ namespace storm {
             }
             
             // Finally, compute the transient probabilities.
-            return this->computeTransientProbabilities<true>(uniformizedMatrix, timeBound, uniformizationRate, totalRewardVector, *this->linearEquationSolver);
+            return this->computeTransientProbabilities<true>(uniformizedMatrix, timeBound, uniformizationRate, totalRewardVector, *this->linearEquationSolverFactory);
         }
         
         template<class ValueType>
@@ -415,7 +416,7 @@ namespace storm {
                 }
             }
             
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseDtmcPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(probabilityMatrix, modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *linearEquationSolver, qualitative)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseDtmcPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(probabilityMatrix, modifiedStateRewardVector, this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *linearEquationSolverFactory, qualitative)));
         }
         
         // Explicitly instantiate the model checker.
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
index 2257075e7..e604e04d4 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h
@@ -3,6 +3,7 @@
 
 #include "src/modelchecker/propositional/SparsePropositionalModelChecker.h"
 #include "src/models/sparse/Ctmc.h"
+#include "src/utility/solver.h"
 #include "src/solver/LinearEquationSolver.h"
 
 namespace storm {
@@ -12,7 +13,7 @@ namespace storm {
         class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker<ValueType> {
         public:
             explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model);
-            explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>&& linearEquationSolver);
+            explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc<ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
             
             // The implemented methods of the AbstractModelChecker interface.
             virtual bool canHandle(storm::logic::Formula const& formula) const override;
@@ -29,7 +30,7 @@ namespace storm {
         private:
             // The methods that perform the actual checking.
             std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::vector<ValueType> const& exitRates, bool qualitative, double lowerBound, double upperBound) const;
-            static std::vector<ValueType> computeUntilProbabilitiesHelper(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::LinearEquationSolver<ValueType> const& linearEquationSolver);
+            static std::vector<ValueType> computeUntilProbabilitiesHelper(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> computeInstantaneousRewardsHelper(double timeBound) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(double timeBound) const;
             std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::BitVector const& targetStates, bool qualitative) const;
@@ -53,13 +54,13 @@ namespace storm {
              * @param timeBound The time bound to use.
              * @param uniformizationRate The used uniformization rate.
              * @param values A vector mapping each state to an initial probability.
-             * @param linearEquationSolver The linear equation solver to use.
+             * @param linearEquationSolverFactory The factory to use when instantiating new linear equation solvers.
              * @param useMixedPoissonProbabilities If set to true, instead of taking the poisson probabilities,  mixed
              * poisson probabilities are used.
              * @return The vector of transient probabilities.
              */
             template<bool useMixedPoissonProbabilities = false>
-            std::vector<ValueType> computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver) const;
+            std::vector<ValueType> computeTransientProbabilities(storm::storage::SparseMatrix<ValueType> const& uniformizedMatrix, ValueType timeBound, ValueType uniformizationRate, std::vector<ValueType> values, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory) const;
             
             /*!
              * Converts the given rate-matrix into a time-abstract probability matrix.
@@ -70,7 +71,7 @@ namespace storm {
             static storm::storage::SparseMatrix<ValueType> computeProbabilityMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates);
             
             // An object that is used for solving linear equations and performing matrix-vector multiplication.
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linearEquationSolver;
+            std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
         
     } // namespace modelchecker
diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
index 45558db78..45dd66b00 100644
--- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
@@ -21,12 +21,12 @@
 namespace storm {
     namespace modelchecker {
         template<typename ValueType>
-        SparseMarkovAutomatonCslModelChecker<ValueType>::SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver) : model(model), nondeterministicLinearEquationSolver(nondeterministicLinearEquationSolver) {
+        SparseMarkovAutomatonCslModelChecker<ValueType>::SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType>>&& nondeterministicLinearEquationSolverFactory) : model(model), nondeterministicLinearEquationSolverFactory(std::move(nondeterministicLinearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
-        SparseMarkovAutomatonCslModelChecker<ValueType>::SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model) : model(model), nondeterministicLinearEquationSolver(storm::utility::solver::getNondeterministicLinearEquationSolver<ValueType>()) {
+        SparseMarkovAutomatonCslModelChecker<ValueType>::SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model) : model(model), nondeterministicLinearEquationSolverFactory(new storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType>()) {
             // Intentionally left empty.
         }
         
@@ -36,7 +36,7 @@ namespace storm {
         }
 
         template<typename ValueType>
-        void SparseMarkovAutomatonCslModelChecker<ValueType>::computeBoundedReachabilityProbabilities(bool min, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps) {
+        void SparseMarkovAutomatonCslModelChecker<ValueType>::computeBoundedReachabilityProbabilities(bool min, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType> const& nondeterministicLinearEquationSolverFactory) {
             // Start by computing four sparse matrices:
             // * a matrix aMarkovian with all (discretized) transitions from Markovian non-goal states to all Markovian non-goal states.
             // * a matrix aMarkovianToProbabilistic with all (discretized) transitions from Markovian non-goal states to all probabilistic non-goal states.
@@ -89,8 +89,8 @@ namespace storm {
                 }
             }
             
-            std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministiclinearEquationSolver = storm::utility::solver::getNondeterministicLinearEquationSolver<ValueType>();
-            
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory.create(aProbabilistic);
+                        
             // Perform the actual value iteration
             // * loop until the step bound has been reached
             // * in the loop:
@@ -106,7 +106,7 @@ namespace storm {
                 storm::utility::vector::addVectorsInPlace(bProbabilistic, bProbabilisticFixed);
                 
                 // Now perform the inner value iteration for probabilistic states.
-                nondeterministiclinearEquationSolver->solveEquationSystem(min, aProbabilistic, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
+                solver->solveEquationSystem(min, aProbabilistic, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
                 
                 // (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic.
                 aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian);
@@ -119,7 +119,7 @@ namespace storm {
             // After the loop, perform one more step of the value iteration for PS states.
             aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
             storm::utility::vector::addVectorsInPlace(bProbabilistic, bProbabilisticFixed);
-            nondeterministiclinearEquationSolver->solveEquationSystem(min, aProbabilistic, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
+            solver->solveEquationSystem(min, aProbabilistic, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
         }
         
         template<typename ValueType>
@@ -149,7 +149,7 @@ namespace storm {
             std::vector<ValueType> vProbabilistic(probabilisticNonGoalStates.getNumberOfSetBits());
             std::vector<ValueType> vMarkovian(markovianNonGoalStates.getNumberOfSetBits());
             
-            computeBoundedReachabilityProbabilities(minimize, transitionMatrix, exitRates, markovianStates, psiStates, markovianNonGoalStates, probabilisticNonGoalStates, vMarkovian, vProbabilistic, delta, numberOfSteps);
+            computeBoundedReachabilityProbabilities(minimize, transitionMatrix, exitRates, markovianStates, psiStates, markovianNonGoalStates, probabilisticNonGoalStates, vMarkovian, vProbabilistic, delta, numberOfSteps, *nondeterministicLinearEquationSolverFactory);
             
             // (4) If the lower bound of interval was non-zero, we need to take the current values as the starting values for a subsequent value iteration.
             if (lowerBound != storm::utility::zero<ValueType>()) {
@@ -167,7 +167,7 @@ namespace storm {
                 STORM_LOG_INFO("Performing " << numberOfSteps << " iterations (delta=" << delta << ") for interval [0, " << lowerBound << "]." << std::endl);
                 
                 // Compute the bounded reachability for interval [0, b-a].
-                computeBoundedReachabilityProbabilities(minimize, transitionMatrix, exitRates, markovianStates, storm::storage::BitVector(model.getNumberOfStates()), markovianStates, ~markovianStates, vAllMarkovian, vAllProbabilistic, delta, numberOfSteps);
+                computeBoundedReachabilityProbabilities(minimize, transitionMatrix, exitRates, markovianStates, storm::storage::BitVector(model.getNumberOfStates()), markovianStates, ~markovianStates, vAllMarkovian, vAllProbabilistic, delta, numberOfSteps, *nondeterministicLinearEquationSolverFactory);
                 
                 // Create the result vector out of vAllProbabilistic and vAllMarkovian and return it.
                 std::vector<ValueType> result(model.getNumberOfStates());
@@ -198,7 +198,7 @@ namespace storm {
         
         template<typename ValueType>
         std::vector<ValueType> SparseMarkovAutomatonCslModelChecker<ValueType>::computeUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const {
-            return storm::modelchecker::SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(minimize, model.getTransitionMatrix(), model.getBackwardTransitions(), phiStates, psiStates, nondeterministicLinearEquationSolver, qualitative);
+            return storm::modelchecker::SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(minimize, model.getTransitionMatrix(), model.getBackwardTransitions(), phiStates, psiStates, *nondeterministicLinearEquationSolverFactory, qualitative);
         }
         
         template<typename ValueType>
@@ -379,7 +379,8 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> sspMatrix = sspMatrixBuilder.build(currentChoice);
             
             std::vector<ValueType> x(numberOfStatesNotInMecs + mecDecomposition.size());
-            nondeterministicLinearEquationSolver->solveEquationSystem(minimize, sspMatrix, x, b);
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory->create(sspMatrix);
+            solver->solveEquationSystem(minimize, sspMatrix, x, b);
             
             // Prepare result vector.
             std::vector<ValueType> result(model.getNumberOfStates());
@@ -437,7 +438,8 @@ namespace storm {
                 
         template<typename ValueType>
         ValueType SparseMarkovAutomatonCslModelChecker<ValueType>::computeLraForMaximalEndComponent(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::BitVector const& markovianStates, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& psiStates, storm::storage::MaximalEndComponent const& mec, uint_fast64_t mecIndex) {
-            std::shared_ptr<storm::solver::LpSolver> solver = storm::utility::solver::getLpSolver("LRA for MEC");
+            std::unique_ptr<storm::utility::solver::LpSolverFactory> lpSolverFactory(new storm::utility::solver::LpSolverFactory());
+            std::unique_ptr<storm::solver::LpSolver> solver = lpSolverFactory->create("LRA for MEC");
             solver->setModelSense(minimize ? storm::solver::LpSolver::ModelSense::Maximize : storm::solver::LpSolver::ModelSense::Minimize);
             
             // First, we need to create the variables for the problem.
@@ -561,8 +563,8 @@ namespace storm {
             storm::utility::vector::selectVectorValuesRepeatedly(b, maybeStates, model.getNondeterministicChoiceIndices(), rewardValues);
             
             // Solve the corresponding system of equations.
-            std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministiclinearEquationSolver = storm::utility::solver::getNondeterministicLinearEquationSolver<ValueType>();
-            nondeterministiclinearEquationSolver->solveEquationSystem(minimize, submatrix, x, b);
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory->create(submatrix);
+            solver->solveEquationSystem(minimize, submatrix, x, b);
             
             // Create resulting vector.
             std::vector<ValueType> result(model.getNumberOfStates());
diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
index c198be08a..4c29dc743 100644
--- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
+++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
@@ -14,7 +14,7 @@ namespace storm {
         template<typename ValueType>
         class SparseMarkovAutomatonCslModelChecker : public AbstractModelChecker {
         public:
-            explicit SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver);
+            explicit SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType>>&& nondeterministicLinearEquationSolver);
             explicit SparseMarkovAutomatonCslModelChecker(storm::models::sparse::MarkovAutomaton<ValueType> const& model);
             
             // The implemented methods of the AbstractModelChecker interface.
@@ -30,7 +30,7 @@ namespace storm {
         private:
             // The methods that perform the actual checking.
             std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& psiStates, std::pair<double, double> const& boundsPair) const;
-            static void computeBoundedReachabilityProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps);
+            static void computeBoundedReachabilityProbabilities(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, std::vector<ValueType> const& exitRates, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& goalStates, storm::storage::BitVector const& markovianNonGoalStates, storm::storage::BitVector const& probabilisticNonGoalStates, std::vector<ValueType>& markovianNonGoalValues, std::vector<ValueType>& probabilisticNonGoalValues, ValueType delta, uint_fast64_t numberOfSteps, storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType> const& nondeterministicLinearEquationSolverFactory);
             std::vector<ValueType> computeUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const;
             std::vector<ValueType> computeReachabilityRewardsHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const;
             std::vector<ValueType> computeLongRunAverageHelper(bool minimize, storm::storage::BitVector const& psiStates, bool qualitative) const;
@@ -66,8 +66,8 @@ namespace storm {
             // The model this model checker is supposed to analyze.
             storm::models::sparse::MarkovAutomaton<ValueType> const& model;
             
-            // A solver that is used for solving systems of linear equations that are the result of nondeterministic choices.
-            std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver;
+            // An object that is used for retrieving solvers for systems of linear equations that are the result of nondeterministic choices.
+            std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType>> nondeterministicLinearEquationSolverFactory;
         };
     }
 }
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index cdde1e716..1b543b5fc 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -16,12 +16,12 @@
 namespace storm {
     namespace modelchecker {
         template<typename ValueType>
-        SparseDtmcPrctlModelChecker<ValueType>::SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>&& linearEquationSolver) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolver(std::move(linearEquationSolver)) {
+        SparseDtmcPrctlModelChecker<ValueType>::SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolverFactory(std::move(linearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
-        SparseDtmcPrctlModelChecker<ValueType>::SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolver(storm::utility::solver::getLinearEquationSolver<ValueType>()) {
+        SparseDtmcPrctlModelChecker<ValueType>::SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model) : SparsePropositionalModelChecker<ValueType>(model), linearEquationSolverFactory(new storm::utility::solver::LinearEquationSolverFactory<ValueType>()) {
             // Intentionally left empty.
         }
         
@@ -53,8 +53,9 @@ namespace storm {
                 storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::one<ValueType>());
                 
                 // Perform the matrix vector multiplication as often as required by the formula bound.
-                STORM_LOG_THROW(linearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-                this->linearEquationSolver->performMatrixVectorMultiplication(submatrix, subresult, nullptr, stepBound);
+                STORM_LOG_THROW(linearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(submatrix);
+                solver->performMatrixVectorMultiplication(subresult, nullptr, stepBound);
                 
                 // Set the values of the resulting vector accordingly.
                 storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, subresult);
@@ -76,13 +77,14 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver) {
+        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::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>());
             
             // Perform one single matrix-vector multiplication.
-            linearEquationSolver.performMatrixVectorMultiplication(transitionMatrix, result);
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix);
+            solver->performMatrixVectorMultiplication(result);
             return result;
         }
         
@@ -90,11 +92,11 @@ namespace storm {
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<ValueType>::computeNextProbabilities(storm::logic::NextFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(pathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeNextProbabilitiesHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolver)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeNextProbabilitiesHelper(this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory)));
         }
         
         template<typename ValueType>
-        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(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::LinearEquationSolver<ValueType> const& linearEquationSolver) {
+        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(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) {
             // 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);
@@ -135,7 +137,8 @@ namespace storm {
                     std::vector<ValueType> b = transitionMatrix.getConstrainedRowSumVector(maybeStates, statesWithProbability1);
                     
                     // Now solve the created system of linear equations.
-                    linearEquationSolver.solveEquationSystem(submatrix, x, b);
+                    std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(submatrix);
+                    solver->solveEquationSystem(x, b);
                     
                     // Set values of resulting vector according to result.
                     storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
@@ -155,7 +158,7 @@ namespace storm {
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();;
             ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();;
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolver)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), qualitative, *this->linearEquationSolverFactory)));
         }
         
         template<typename ValueType>
@@ -183,8 +186,9 @@ namespace storm {
             }
             
             // Perform the matrix vector multiplication as often as required by the formula bound.
-            STORM_LOG_THROW(linearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-            this->linearEquationSolver->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), result, &totalRewardVector, stepBound);
+            STORM_LOG_THROW(linearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
+            solver->performMatrixVectorMultiplication(result, &totalRewardVector, stepBound);
             
             return result;
         }
@@ -204,8 +208,9 @@ namespace storm {
             std::vector<ValueType> result(this->getModel().getStateRewardVector());
             
             // Perform the matrix vector multiplication as often as required by the formula bound.
-            STORM_LOG_THROW(linearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-            this->linearEquationSolver->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), result, nullptr, stepCount);
+            STORM_LOG_THROW(linearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
+            solver->performMatrixVectorMultiplication(result, nullptr, stepCount);
             
             return result;
         }
@@ -217,7 +222,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver, bool qualitative) {
+        std::vector<ValueType> SparseDtmcPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative) {
             // Only compute the result if the model has at least one reward this->getModel().
             STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
             
@@ -278,7 +283,8 @@ namespace storm {
                 }
                 
                 // Now solve the resulting equation system.
-                linearEquationSolver.solveEquationSystem(submatrix, x, b);
+                std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(submatrix);
+                solver->solveEquationSystem(x, b);
                 
                 // Set values of resulting vector according to result.
                 storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
@@ -295,7 +301,7 @@ namespace storm {
         std::unique_ptr<CheckResult> SparseDtmcPrctlModelChecker<ValueType>::computeReachabilityRewards(storm::logic::ReachabilityRewardFormula const& rewardPathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->linearEquationSolver, qualitative)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative)));
         }
         
         template<typename ValueType>
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
index b5e73436b..e5407f9b5 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
@@ -17,7 +17,7 @@ namespace storm {
             friend class SparseCtmcCslModelChecker<ValueType>;
             
             explicit SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model);
-            explicit SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model, std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>&& linearEquationSolver);
+            explicit SparseDtmcPrctlModelChecker(storm::models::sparse::Dtmc<ValueType> const& model, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory);
             
             // The implemented methods of the AbstractModelChecker interface.
             virtual bool canHandle(storm::logic::Formula const& formula) const override;
@@ -34,14 +34,14 @@ namespace storm {
         private:
             // The methods that perform the actual checking.
             std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound) const;
-            static std::vector<ValueType> computeNextProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver);
-            static std::vector<ValueType> computeUntilProbabilitiesHelper(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::LinearEquationSolver<ValueType> const& linearEquationSolver);
+            static std::vector<ValueType> computeNextProbabilitiesHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory);
+            static std::vector<ValueType> computeUntilProbabilitiesHelper(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> computeInstantaneousRewardsHelper(uint_fast64_t stepCount) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(uint_fast64_t stepBound) const;
-            static std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::solver::LinearEquationSolver<ValueType> const& linearEquationSolver, bool qualitative);
+            static std::vector<ValueType> computeReachabilityRewardsHelper(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory<ValueType> const& linearEquationSolverFactory, bool qualitative);
             
-            // An object that is used for solving linear equations and performing matrix-vector multiplication.
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> linearEquationSolver;
+            // An object that is used for retrieving linear equation solvers.
+            std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<ValueType>> linearEquationSolverFactory;
         };
         
     } // namespace modelchecker
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index fdf85d1b6..8b5975aca 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -16,12 +16,12 @@
 namespace storm {
     namespace modelchecker {
         template<typename ValueType>
-        SparseMdpPrctlModelChecker<ValueType>::SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model) : SparsePropositionalModelChecker<ValueType>(model), nondeterministicLinearEquationSolver(storm::utility::solver::getNondeterministicLinearEquationSolver<ValueType>()) {
+        SparseMdpPrctlModelChecker<ValueType>::SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model) : SparsePropositionalModelChecker<ValueType>(model), nondeterministicLinearEquationSolverFactory(new storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType>()) {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
-        SparseMdpPrctlModelChecker<ValueType>::SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver) : SparsePropositionalModelChecker<ValueType>(model), nondeterministicLinearEquationSolver(nondeterministicLinearEquationSolver) {
+        SparseMdpPrctlModelChecker<ValueType>::SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType>>&& nondeterministicLinearEquationSolverFactory) : SparsePropositionalModelChecker<ValueType>(model), nondeterministicLinearEquationSolverFactory(std::move(nondeterministicLinearEquationSolverFactory)) {
             // Intentionally left empty.
         }
         
@@ -57,8 +57,9 @@ namespace storm {
                 std::vector<ValueType> subresult(statesWithProbabilityGreater0.getNumberOfSetBits());
                 storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::one<ValueType>());
             
-                STORM_LOG_THROW(nondeterministicLinearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid equation solver available.");
-                this->nondeterministicLinearEquationSolver->performMatrixVectorMultiplication(minimize, submatrix, subresult, nullptr, stepBound);
+                STORM_LOG_THROW(nondeterministicLinearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid equation solver available.");
+                std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory->create(submatrix);
+                solver->performMatrixVectorMultiplication(minimize, subresult, nullptr, stepBound);
                 
                 // Set the values of the resulting vector accordingly.
                 storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, subresult);
@@ -85,8 +86,9 @@ namespace storm {
             std::vector<ValueType> result(this->getModel().getNumberOfStates());
             storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
             
-            STORM_LOG_THROW(nondeterministicLinearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid equation solver available.");
-            this->nondeterministicLinearEquationSolver->performMatrixVectorMultiplication(minimize, this->getModel().getTransitionMatrix(), result);
+            STORM_LOG_THROW(nondeterministicLinearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid equation solver available.");
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
+            solver->performMatrixVectorMultiplication(minimize, result);
             
             return result;
         }
@@ -101,11 +103,11 @@ namespace storm {
         
         template<typename ValueType>
         std::vector<ValueType> SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const {
-            return computeUntilProbabilitiesHelper(minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, nondeterministicLinearEquationSolver, qualitative);
+            return computeUntilProbabilitiesHelper(minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), phiStates, psiStates, *nondeterministicLinearEquationSolverFactory, qualitative);
         }
         
         template<typename ValueType>
-        std::vector<ValueType> SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver, bool qualitative) {
+        std::vector<ValueType> SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType> const& nondeterministicLinearEquationSolverFactory, bool qualitative) {
             size_t numberOfStates = phiStates.size();
             
             // We need to identify the states which have to be taken out of the matrix, i.e.
@@ -146,7 +148,8 @@ namespace storm {
                     std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
                     
                     // Solve the corresponding system of equations.
-                    nondeterministicLinearEquationSolver->solveEquationSystem(minimize, submatrix, x, b);
+                    std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory.create(submatrix);
+                    solver->solveEquationSystem(minimize, x, b);
                     
                     // Set values of resulting vector according to result.
                     storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
@@ -167,7 +170,7 @@ namespace storm {
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();
             ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), nondeterministicLinearEquationSolver, qualitative)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(SparseMdpPrctlModelChecker<ValueType>::computeUntilProbabilitiesHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), *nondeterministicLinearEquationSolverFactory, qualitative)));
         }
         
         template<typename ValueType>
@@ -194,7 +197,8 @@ namespace storm {
                 result.resize(this->getModel().getNumberOfStates());
             }
             
-            this->nondeterministicLinearEquationSolver->performMatrixVectorMultiplication(minimize, this->getModel().getTransitionMatrix(), result, &totalRewardVector, stepBound);
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
+            solver->performMatrixVectorMultiplication(minimize, result, &totalRewardVector, stepBound);
             
             return result;
         }
@@ -214,8 +218,9 @@ namespace storm {
             // Initialize result to state rewards of the this->getModel().
             std::vector<ValueType> result(this->getModel().getStateRewardVector());
             
-            STORM_LOG_THROW(nondeterministicLinearEquationSolver != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
-            this->nondeterministicLinearEquationSolver->performMatrixVectorMultiplication(minimize, this->getModel().getTransitionMatrix(), result, nullptr, stepCount);
+            STORM_LOG_THROW(nondeterministicLinearEquationSolverFactory != nullptr, storm::exceptions::InvalidStateException, "No valid linear equation solver available.");
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory->create(this->getModel().getTransitionMatrix());
+            solver->performMatrixVectorMultiplication(minimize, result, nullptr, stepCount);
             
             return result;
         }
@@ -228,7 +233,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::vector<ValueType> SparseMdpPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::solver::NondeterministicLinearEquationSolver<ValueType> const& nondeterministicLinearEquationSolver, bool qualitative) const {
+        std::vector<ValueType> SparseMdpPrctlModelChecker<ValueType>::computeReachabilityRewardsHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType> const& nondeterministicLinearEquationSolverFactory, bool qualitative) const {
             // Only compute the result if the model has at least one reward this->getModel().
             STORM_LOG_THROW(stateRewardVector || transitionRewardMatrix, storm::exceptions::InvalidPropertyException, "Missing reward model for formula. Skipping formula.");
             
@@ -294,7 +299,8 @@ namespace storm {
                 std::vector<ValueType> x(maybeStates.getNumberOfSetBits());
                 
                 // Solve the corresponding system of equations.
-                this->nondeterministicLinearEquationSolver->solveEquationSystem(minimize, submatrix, x, b);
+                std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> solver = nondeterministicLinearEquationSolverFactory.create(submatrix);
+                solver->solveEquationSystem(minimize, x, b);
                 
                 // Set values of resulting vector according to result.
                 storm::utility::vector::setVectorValues<ValueType>(result, maybeStates, x);
@@ -312,7 +318,7 @@ namespace storm {
             STORM_LOG_THROW(optimalityType, storm::exceptions::InvalidArgumentException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic.");
             std::unique_ptr<CheckResult> subResultPointer = this->check(rewardPathFormula.getSubformula());
             ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult();
-            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->nondeterministicLinearEquationSolver, qualitative)));
+            return std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeReachabilityRewardsHelper(optimalityType.get() == storm::logic::OptimalityType::Minimize, this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->nondeterministicLinearEquationSolverFactory, qualitative)));
         }
         
         template<typename ValueType>
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
index 153b5168a..b712f1f2f 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
@@ -29,7 +29,7 @@ namespace storm {
             friend class storm::counterexamples::MILPMinimalLabelSetGenerator<ValueType>;
             
             explicit SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model);
-            explicit SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver);
+            explicit SparseMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType>>&& nondeterministicLinearEquationSolverFactory);
             
             // The implemented methods of the AbstractModelChecker interface.
             virtual bool canHandle(storm::logic::Formula const& formula) const override;
@@ -48,13 +48,13 @@ namespace storm {
             std::vector<ValueType> computeBoundedUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound) const;
             std::vector<ValueType> computeNextProbabilitiesHelper(bool minimize, storm::storage::BitVector const& nextStates);
             std::vector<ValueType> computeUntilProbabilitiesHelper(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const;
-            static std::vector<ValueType> computeUntilProbabilitiesHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver, bool qualitative);
+            static std::vector<ValueType> computeUntilProbabilitiesHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType> const& nondeterministicLinearEquationSolverFactory, bool qualitative);
             std::vector<ValueType> computeInstantaneousRewardsHelper(bool minimize, uint_fast64_t stepCount) const;
             std::vector<ValueType> computeCumulativeRewardsHelper(bool minimize, uint_fast64_t stepBound) const;
-            std::vector<ValueType> computeReachabilityRewardsHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::solver::NondeterministicLinearEquationSolver<ValueType> const& nondeterministicLinearEquationSolver, bool qualitative) const;
+            std::vector<ValueType> computeReachabilityRewardsHelper(bool minimize, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, boost::optional<std::vector<ValueType>> const& stateRewardVector, boost::optional<storm::storage::SparseMatrix<ValueType>> const& transitionRewardMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType> const& nondeterministicLinearEquationSolverFactory, bool qualitative) const;
             
-            // A solver that is used for solving systems of linear equations that are the result of nondeterministic choices.
-            std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver;
+            // An object that is used for retrieving solvers for systems of linear equations that are the result of nondeterministic choices.
+            std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<ValueType>> nondeterministicLinearEquationSolverFactory;
         };
     } // namespace modelchecker
 } // namespace storm
diff --git a/src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h b/src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h
deleted file mode 100644
index 62fd98bb0..000000000
--- a/src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h
+++ /dev/null
@@ -1,41 +0,0 @@
-#ifndef STORM_MODELCHECKER_PRCTL_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_
-#define STORM_MODELCHECKER_PRCTL_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_
-
-#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
-#include "src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.h"
-#include "src/exceptions/InvalidPropertyException.h"
-#include <cmath>
-
-namespace storm {
-    namespace modelchecker {
-        
-        /*
-         * An implementation of the SparseMdpPrctlModelChecker interface that uses topoligical value iteration for solving
-         * equation systems.
-         */
-        template <class ValueType>
-        class TopologicalValueIterationMdpPrctlModelChecker : public SparseMdpPrctlModelChecker<ValueType> {
-        public:
-            /*!
-             * Constructs a SparseMdpPrctlModelChecker with the given model.
-             *
-             * @param model The MDP to be checked.
-             */
-            explicit TopologicalValueIterationMdpPrctlModelChecker(storm::models::sparse::Mdp<ValueType> const& model) : SparseMdpPrctlModelChecker<ValueType>(model, std::shared_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>>(new storm::solver::TopologicalValueIterationNondeterministicLinearEquationSolver<ValueType>())) {
-                // Intentionally left empty.
-            }
-            
-            /*!
-             * Copy constructs a SparseMdpPrctlModelChecker from the given model checker. In particular, this means that the newly
-             * constructed model checker will have the model of the given model checker as its associated model.
-             */
-            explicit TopologicalValueIterationMdpPrctlModelChecker(storm::modelchecker::TopologicalValueIterationMdpPrctlModelChecker<ValueType> const& modelchecker)
-            : SparseMdpPrctlModelChecker<ValueType>(modelchecker) {
-                // Intentionally left empty.
-            }
-        };
-        
-    } // namespace modelchecker
-} // namespace storm
-
-#endif /* STORM_MODELCHECKER_PRCTL_TOPOLOGICALVALUEITERATIONSMDPPRCTLMODELCHECKER_H_ */
diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp
index 63bb8f224..d40e9bb70 100644
--- a/src/solver/GmmxxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxLinearEquationSolver.cpp
@@ -9,19 +9,18 @@
 #include "src/utility/constants.h"
 #include "src/exceptions/InvalidStateException.h"
 
-#include "gmm/gmm_matrix.h"
 #include "gmm/gmm_iter_solvers.h"
 
 namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, bool relative, uint_fast64_t restart) : method(method), precision(precision), maximalNumberOfIterations(maximalNumberOfIterations), preconditioner(preconditioner), relative(relative), restart(restart) {
+        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() {
+        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : originalA(&A), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)) {
             // Get the settings object to customize linear solving.
             storm::settings::modules::GmmxxEquationSolverSettings const& settings = storm::settings::gmmxxEquationSolverSettings();
             
@@ -55,46 +54,39 @@ namespace storm {
         }
         
         template<typename ValueType>
-        LinearEquationSolver<ValueType>* GmmxxLinearEquationSolver<ValueType>::clone() const {
-            return new GmmxxLinearEquationSolver<ValueType>(*this);
-        }
-        
-        template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+        void GmmxxLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             LOG4CPLUS_INFO(logger, "Using method '" << methodToString() << "' with preconditioner " << preconditionerToString() << "'.");
             if (method == SolutionMethod::Jacobi && preconditioner != Preconditioner::None) {
                 LOG4CPLUS_WARN(logger, "Jacobi method currently does not support preconditioners. The requested preconditioner will be ignored.");
             }
             
             if (method == SolutionMethod::Bicgstab || method == SolutionMethod::Qmr || method == SolutionMethod::Gmres) {
-                std::unique_ptr<gmm::csr_matrix<ValueType>> gmmA = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
-                
                 // Prepare an iteration object that determines the accuracy and the maximum number of iterations.
                 gmm::iteration iter(precision, 0, maximalNumberOfIterations);
                 
                 if (method == SolutionMethod::Bicgstab) {
                     if (preconditioner == Preconditioner::Ilu) {
-                        gmm::bicgstab(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmA), iter);
+                        gmm::bicgstab(*gmmxxMatrix, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), iter);
                     } else if (preconditioner == Preconditioner::Diagonal) {
-                        gmm::bicgstab(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmA), iter);
+                        gmm::bicgstab(*gmmxxMatrix, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), iter);
                     } else if (preconditioner == Preconditioner::None) {
-                        gmm::bicgstab(*gmmA, x, b, gmm::identity_matrix(), iter);
+                        gmm::bicgstab(*gmmxxMatrix, x, b, gmm::identity_matrix(), iter);
                     }
                 } else if (method == SolutionMethod::Qmr) {
                     if (preconditioner == Preconditioner::Ilu) {
-                        gmm::qmr(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmA), iter);
+                        gmm::qmr(*gmmxxMatrix, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), iter);
                     } else if (preconditioner == Preconditioner::Diagonal) {
-                        gmm::qmr(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmA), iter);
+                        gmm::qmr(*gmmxxMatrix, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), iter);
                     } else if (preconditioner == Preconditioner::None) {
-                        gmm::qmr(*gmmA, x, b, gmm::identity_matrix(), iter);
+                        gmm::qmr(*gmmxxMatrix, x, b, gmm::identity_matrix(), iter);
                     }
                 } else if (method == SolutionMethod::Gmres) {
                     if (preconditioner == Preconditioner::Ilu) {
-                        gmm::gmres(*gmmA, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmA), restart, iter);
+                        gmm::gmres(*gmmxxMatrix, x, b, gmm::ilu_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), restart, iter);
                     } else if (preconditioner == Preconditioner::Diagonal) {
-                        gmm::gmres(*gmmA, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmA), restart, iter);
+                        gmm::gmres(*gmmxxMatrix, x, b, gmm::diagonal_precond<gmm::csr_matrix<ValueType>>(*gmmxxMatrix), restart, iter);
                     } else if (preconditioner == Preconditioner::None) {
-                        gmm::gmres(*gmmA, x, b, gmm::identity_matrix(), restart, iter);
+                        gmm::gmres(*gmmxxMatrix, x, b, gmm::identity_matrix(), restart, iter);
                     }
                 }
                 
@@ -105,7 +97,7 @@ namespace storm {
                     LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
                 }
             } else if (method == SolutionMethod::Jacobi) {
-                uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b, multiplyResult);
+                uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(*originalA, x, b, multiplyResult);
                 
                 // Check if the solver converged and issue a warning otherwise.
                 if (iterations < maximalNumberOfIterations) {
@@ -117,10 +109,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
-            // Transform the transition probability A to the gmm++ format to use its arithmetic.
-            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
-
+        void GmmxxLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
             // Set up some temporary variables so that we can just swap pointers instead of copying the result after
             // each iteration.
             std::vector<ValueType>* currentX = &x;
diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h
index ac3adfdb2..7ec0dcc8f 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 "gmm/gmm_matrix.h"
+
 #include "LinearEquationSolver.h"
 
 namespace storm {
@@ -25,6 +27,7 @@ namespace storm {
             /*!
              * 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.
@@ -34,18 +37,18 @@ namespace storm {
              * @param restart An optional argument that specifies after how many iterations restarted methods are
              * supposed to actually to a restart.
              */
-            GmmxxLinearEquationSolver(SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, Preconditioner preconditioner, bool relative = true, uint_fast64_t restart = 0);
+            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);
             
             /*!
              * 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();
-            
-            virtual LinearEquationSolver<ValueType>* clone() const override;
+            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 solveEquationSystem(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            
-            virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
         private:
             /*!
@@ -76,6 +79,12 @@ namespace storm {
              */
             std::string preconditionerToString() const;
             
+            // A pointer to the original sparse matrix given to this solver.
+            storm::storage::SparseMatrix<ValueType> const* originalA;
+            
+            // 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;
 
diff --git a/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp b/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp
index abd11ad73..94fad460c 100644
--- a/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp
+++ b/src/solver/GmmxxNondeterministicLinearEquationSolver.cpp
@@ -10,7 +10,7 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        GmmxxNondeterministicLinearEquationSolver<ValueType>::GmmxxNondeterministicLinearEquationSolver() {
+        GmmxxNondeterministicLinearEquationSolver<ValueType>::GmmxxNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), rowGroupIndices(A.getRowGroupIndices()) {
             // Get the settings object to customize solving.
             storm::settings::modules::GmmxxEquationSolverSettings const& settings = storm::settings::gmmxxEquationSolverSettings();
             
@@ -21,25 +21,17 @@ namespace storm {
         }
         
         template<typename ValueType>
-        GmmxxNondeterministicLinearEquationSolver<ValueType>::GmmxxNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
+        GmmxxNondeterministicLinearEquationSolver<ValueType>::GmmxxNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), rowGroupIndices(A.getRowGroupIndices()), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
             // Intentionally left empty.
         }
 
         
         template<typename ValueType>
-        NondeterministicLinearEquationSolver<ValueType>* GmmxxNondeterministicLinearEquationSolver<ValueType>::clone() const {
-            return new GmmxxNondeterministicLinearEquationSolver<ValueType>(*this);
-        }
-        
-        template<typename ValueType>
-        void GmmxxNondeterministicLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
-            // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
-            
+        void GmmxxNondeterministicLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
             // Set up the environment for the power method. If scratch memory was not provided, we need to create it.
             bool multiplyResultMemoryProvided = true;
             if (multiplyResult == nullptr) {
-                multiplyResult = new std::vector<ValueType>(A.getRowCount());
+                multiplyResult = new std::vector<ValueType>(b.size());
                 multiplyResultMemoryProvided = false;
             }
             
@@ -64,9 +56,9 @@ namespace storm {
                 
                 // Reduce the vector x by applying min/max over all nondeterministic choices.
                 if (minimize) {
-                    storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, A.getRowGroupIndices());
+                    storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, rowGroupIndices);
                 } else {
-                    storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, A.getRowGroupIndices());
+                    storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, rowGroupIndices);
                 }
                 
                 // Determine whether the method converged.
@@ -100,13 +92,10 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void GmmxxNondeterministicLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
-            // Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
-            
+        void GmmxxNondeterministicLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(bool minimize, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
             bool multiplyResultMemoryProvided = true;
             if (multiplyResult == nullptr) {
-                multiplyResult = new std::vector<ValueType>(A.getRowCount());
+                multiplyResult = new std::vector<ValueType>(rowGroupIndices.size() - 1);
                 multiplyResultMemoryProvided = false;
             }
             
@@ -119,9 +108,9 @@ namespace storm {
                 }
                 
                 if (minimize) {
-                    storm::utility::vector::reduceVectorMin(*multiplyResult, x, A.getRowGroupIndices());
+                    storm::utility::vector::reduceVectorMin(*multiplyResult, x, rowGroupIndices);
                 } else {
-                    storm::utility::vector::reduceVectorMax(*multiplyResult, x, A.getRowGroupIndices());
+                    storm::utility::vector::reduceVectorMax(*multiplyResult, x, rowGroupIndices);
                 }
             }
             
diff --git a/src/solver/GmmxxNondeterministicLinearEquationSolver.h b/src/solver/GmmxxNondeterministicLinearEquationSolver.h
index 799853173..2d87c60ca 100644
--- a/src/solver/GmmxxNondeterministicLinearEquationSolver.h
+++ b/src/solver/GmmxxNondeterministicLinearEquationSolver.h
@@ -1,6 +1,8 @@
 #ifndef STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_
 #define STORM_SOLVER_GMMXXNONDETERMINISTICLINEAREQUATIONSOLVER_H_
 
+#include "gmm/gmm_matrix.h"
+
 #include "src/solver/NondeterministicLinearEquationSolver.h"
 
 namespace storm {
@@ -9,32 +11,39 @@ namespace storm {
         /*!
          * A class that uses the gmm++ library to implement the NondeterministicLinearEquationSolver interface.
          */
-        template<class Type>
-        class GmmxxNondeterministicLinearEquationSolver : public NondeterministicLinearEquationSolver<Type> {
+        template<class ValueType>
+        class GmmxxNondeterministicLinearEquationSolver : public NondeterministicLinearEquationSolver<ValueType> {
         public:
             /*!
              * Constructs a nondeterministic linear equation solver with parameters being set according to the settings
              * object.
+             *
+             * @param A The matrix defining the coefficients of the linear equation system.             
              */
-            GmmxxNondeterministicLinearEquationSolver();
+            GmmxxNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
             
             /*!
              * Constructs a nondeterminstic linear equation solver with the given parameters.
              *
+             * @param A The matrix defining the coefficients of the linear equation system.
              * @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.
              */
-            GmmxxNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
+            GmmxxNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
 
-            virtual NondeterministicLinearEquationSolver<Type>* clone() const;
-            
-            virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1, std::vector<Type>* multiplyResult = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(bool minimize, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
             
-            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<Type>* multiplyResult = nullptr, std::vector<Type>* newX = nullptr) const override;
+            virtual void solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
             
         private:
+            // The (gmm++) matrix associated with this equation solver.
+            std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix;
+            
+            // A reference to the row group indices of the original matrix.
+            std::vector<uint_fast64_t> const& rowGroupIndices;
+            
             // The required precision for the iterative methods.
             double precision;
             
diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h
index 2a0f6579e..82a229659 100644
--- a/src/solver/LinearEquationSolver.h
+++ b/src/solver/LinearEquationSolver.h
@@ -15,30 +15,23 @@ namespace storm {
         template<class Type>
         class LinearEquationSolver {
         public:
-            /*!
-             * Makes a copy of the linear equation solver.
-             *
-             * @return A pointer to a copy of the linear equation solver.
-             */
-            virtual LinearEquationSolver<Type>* clone() const = 0;
-            
             /*!
              * Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution.
-             * The solution of the set of linear equations will be written to the vector x.
+             * The solution of the set of linear equations will be written to the vector x. Note that the matrix A has
+             * to be given upon construction time of the solver object.
              *
-             * @param A The coefficient matrix of the equation system.
              * @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A.
              * @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A.
              * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
              * vector must be equal to the number of rows of A.
              */
-            virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type> const& b, std::vector<Type>* multiplyResult = nullptr) const = 0;
+            virtual void solveEquationSystem(std::vector<Type>& x, std::vector<Type> const& b, std::vector<Type>* multiplyResult = nullptr) const = 0;
             
             /*!
              * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After
-             * performing the necessary multiplications, the result is written to the input vector x.
+             * performing the necessary multiplications, the result is written to the input vector x. Note that the
+             * matrix A has to be given upon construction time of the solver object.
              *
-             * @param A The matrix to use for matrix-vector multiplication.
              * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal
              * to the number of rows of A.
              * @param b If non-null, this vector is added after each multiplication. If given, its length must be equal
@@ -46,7 +39,7 @@ namespace storm {
              * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
              * vector must be equal to the number of rows of A.
              */
-            virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& A, std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1, std::vector<Type>* multiplyResult = nullptr) const = 0;
+            virtual void performMatrixVectorMultiplication(std::vector<Type>& x, std::vector<Type>* b = nullptr, uint_fast64_t n = 1, std::vector<Type>* multiplyResult = nullptr) const = 0;
         };
         
     } // namespace solver
diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp
index a872643b4..6b89d8134 100644
--- a/src/solver/NativeLinearEquationSolver.cpp
+++ b/src/solver/NativeLinearEquationSolver.cpp
@@ -10,12 +10,12 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : method(method), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod 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() {
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) {
             // Get the settings object to customize linear solving.
             storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
             
@@ -30,14 +30,9 @@ namespace storm {
                 method = SolutionMethod::Jacobi;
             }
         }
-        
-        template<typename ValueType>
-        LinearEquationSolver<ValueType>* NativeLinearEquationSolver<ValueType>::clone() const {
-            return new NativeLinearEquationSolver<ValueType>(*this);
-        }
-        
+                
         template<typename ValueType>
-        void NativeLinearEquationSolver<ValueType>::solveEquationSystem(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+        void NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
             // Get a Jacobi decomposition of the matrix A.
             std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
             
@@ -88,7 +83,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void NativeLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+        void NativeLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
             // Set up some temporary variables so that we can just swap pointers instead of copying the result after
             // each iteration.
             std::vector<ValueType>* currentX = &x;
diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h
index 69504b627..ae008a5aa 100644
--- a/src/solver/NativeLinearEquationSolver.h
+++ b/src/solver/NativeLinearEquationSolver.h
@@ -19,25 +19,26 @@ namespace storm {
             
             /*!
              * 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.
              */
-            NativeLinearEquationSolver();
+            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
 
             /*!
              * 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(SolutionMethod method, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
-            
-            virtual LinearEquationSolver<ValueType>* clone() const override;
+            NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, SolutionMethod 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 solveEquationSystem(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            
-            virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
 
         private:
             /*!
@@ -47,6 +48,9 @@ namespace storm {
              */
             std::string methodToString() const;
             
+            // A reference to the matrix the gives the coefficients of the linear equation system.
+            storm::storage::SparseMatrix<ValueType> const& A;
+            
             // The method to use for solving linear equation systems.
             SolutionMethod method;
             
diff --git a/src/solver/NativeNondeterministicLinearEquationSolver.cpp b/src/solver/NativeNondeterministicLinearEquationSolver.cpp
index 51698aed7..7965c9dda 100644
--- a/src/solver/NativeNondeterministicLinearEquationSolver.cpp
+++ b/src/solver/NativeNondeterministicLinearEquationSolver.cpp
@@ -9,7 +9,7 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        NativeNondeterministicLinearEquationSolver<ValueType>::NativeNondeterministicLinearEquationSolver() {
+        NativeNondeterministicLinearEquationSolver<ValueType>::NativeNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) {
             // Get the settings object to customize solving.
             storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
             
@@ -20,17 +20,12 @@ namespace storm {
         }
         
         template<typename ValueType>
-        NativeNondeterministicLinearEquationSolver<ValueType>::NativeNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
+        NativeNondeterministicLinearEquationSolver<ValueType>::NativeNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : A(A), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations) {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
-        NondeterministicLinearEquationSolver<ValueType>* NativeNondeterministicLinearEquationSolver<ValueType>::clone() const {
-            return new NativeNondeterministicLinearEquationSolver<ValueType>(*this);
-        }
-        
-        template<typename ValueType>
-        void NativeNondeterministicLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+        void NativeNondeterministicLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
             
             // Set up the environment for the power method. If scratch memory was not provided, we need to create it.
             bool multiplyResultMemoryProvided = true;
@@ -96,7 +91,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void NativeNondeterministicLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+        void NativeNondeterministicLinearEquationSolver<ValueType>::performMatrixVectorMultiplication(bool minimize, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
             
             // If scratch memory was not provided, we need to create it.
             bool multiplyResultMemoryProvided = true;
diff --git a/src/solver/NativeNondeterministicLinearEquationSolver.h b/src/solver/NativeNondeterministicLinearEquationSolver.h
index 60dd09149..b70a3b03d 100644
--- a/src/solver/NativeNondeterministicLinearEquationSolver.h
+++ b/src/solver/NativeNondeterministicLinearEquationSolver.h
@@ -15,26 +15,30 @@ namespace storm {
             /*!
              * Constructs a nondeterministic linear equation solver with parameters being set according to the settings
              * object.
+             *
+             * @param A The matrix defining the coefficients of the linear equation system.
              */
-            NativeNondeterministicLinearEquationSolver();
+            NativeNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
             
             /*!
              * Constructs a nondeterminstic linear equation solver with the given parameters.
              *
+             * @param A The matrix defining the coefficients of the linear equation system.
              * @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.
              */
-            NativeNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
-            
-            virtual NondeterministicLinearEquationSolver<ValueType>* clone() const override;
+            NativeNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
             
-            virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* newX = nullptr) const override;
+            virtual void performMatrixVectorMultiplication(bool minimize, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* newX = nullptr) const override;
             
-            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
+            virtual void solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
 
         protected:
+            // A reference to the matrix the gives the coefficients of the linear equation system.
+            storm::storage::SparseMatrix<ValueType> const& A;
+            
             // The required precision for the iterative methods.
             double precision;
             
diff --git a/src/solver/NondeterministicLinearEquationSolver.h b/src/solver/NondeterministicLinearEquationSolver.h
index b90773ee6..28751d87a 100644
--- a/src/solver/NondeterministicLinearEquationSolver.h
+++ b/src/solver/NondeterministicLinearEquationSolver.h
@@ -17,18 +17,11 @@ namespace storm {
         class NondeterministicLinearEquationSolver {
         public:
             /*!
-             * Makes a copy of the nondeterministic linear equation solver.
-             *
-             * @return A pointer to a copy of the nondeterministic linear equation solver.
-             */
-            virtual NondeterministicLinearEquationSolver<ValueType>* clone() const = 0;
-            
-            /*!
-             * Solves the equation system x = min/max(A*x + b) given by the parameters.
+             * Solves the equation system x = min/max(A*x + b) given by the parameters. Note that the matrix A has
+             * to be given upon construction time of the solver object.
              *
              * @param minimize If set, all the value of a group of rows is the taken as the minimum over all rows and as
              * the maximum otherwise.
-             * @param A The matrix specifying the coefficients of the linear equations.
              * @param x The solution vector x. The initial values of x represent a guess of the real values to the
              * solver, but may be ignored.
              * @param b The vector to add after matrix-vector multiplication.
@@ -38,17 +31,17 @@ namespace storm {
              * vector must be equal to the length of the vector x (and thus to the number of columns of A).
              * @return The solution vector x of the system of linear equations as the content of the parameter x.
              */
-            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const = 0;
+            virtual void solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const = 0;
             
             /*!
              * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes
              * x[i+1] = min/max(A*x[i] + b) until x[n], where x[0] = x. After each multiplication and addition, the
              * minimal/maximal value out of each row group is selected to reduce the resulting vector to obtain the
-             * vector for the next iteration.
+             * vector for the next iteration. Note that the matrix A has to be given upon construction time of the
+             * solver object.
              *
              * @param minimize If set, all the value of a group of rows is the taken as the minimum over all rows and as
              * the maximum otherwise.
-             * @param A The matrix that is to be multiplied with the vector.
              * @param x The initial vector that is to be multiplied with the matrix. This is also the output parameter,
              * i.e. after the method returns, this vector will contain the computed values.
              * @param b If not null, this vector is added after each multiplication.
@@ -57,7 +50,7 @@ namespace storm {
              * vector must be equal to the number of rows of A.
              * @return The result of the repeated matrix-vector multiplication as the content of the vector x.
              */
-            virtual void performMatrixVectorMultiplication(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
+            virtual void performMatrixVectorMultiplication(bool minimize, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
         };
         
     } // namespace solver
diff --git a/src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.cpp b/src/solver/TopologicalNondeterministicLinearEquationSolver.cpp
similarity index 88%
rename from src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.cpp
rename to src/solver/TopologicalNondeterministicLinearEquationSolver.cpp
index 0bdf1df80..343f36f1b 100644
--- a/src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.cpp
+++ b/src/solver/TopologicalNondeterministicLinearEquationSolver.cpp
@@ -1,4 +1,4 @@
-#include "src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.h"
+#include "src/solver/TopologicalNondeterministicLinearEquationSolver.h"
 
 #include <utility>
 #include <chrono>
@@ -23,7 +23,7 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-		TopologicalValueIterationNondeterministicLinearEquationSolver<ValueType>::TopologicalValueIterationNondeterministicLinearEquationSolver() {
+        TopologicalNondeterministicLinearEquationSolver<ValueType>::TopologicalNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : NativeNondeterministicLinearEquationSolver<ValueType>(A) {
 			// Get the settings object to customize solving.
 			
 			//storm::settings::Settings* settings = storm::settings::Settings::getInstance();
@@ -44,31 +44,28 @@ namespace storm {
         }
         
         template<typename ValueType>
-		TopologicalValueIterationNondeterministicLinearEquationSolver<ValueType>::TopologicalValueIterationNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : NativeNondeterministicLinearEquationSolver<ValueType>(precision, maximalNumberOfIterations, relative) {
+		TopologicalNondeterministicLinearEquationSolver<ValueType>::TopologicalNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : NativeNondeterministicLinearEquationSolver<ValueType>(A, precision, maximalNumberOfIterations, relative) {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
-		NondeterministicLinearEquationSolver<ValueType>* TopologicalValueIterationNondeterministicLinearEquationSolver<ValueType>::clone() const {
-			return new TopologicalValueIterationNondeterministicLinearEquationSolver<ValueType>(*this);
-        }
-        
-        template<typename ValueType>
-		void TopologicalValueIterationNondeterministicLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+		void TopologicalNondeterministicLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
 			
 #ifdef GPU_USE_FLOAT
 #define __FORCE_FLOAT_CALCULATION true
 #else
 #define __FORCE_FLOAT_CALCULATION false
 #endif
-			if (__FORCE_FLOAT_CALCULATION && (sizeof(ValueType) == sizeof(double))) {
-				TopologicalValueIterationNondeterministicLinearEquationSolver<float> tvindles{ this->precision, this->maximalNumberOfIterations, this->relative };
+            if (__FORCE_FLOAT_CALCULATION && std::is_same<ValueType, double>::value) {
+                // FIXME: This actually allocates quite some storage, because of this conversion, is it really necessary?
+				storm::storage::SparseMatrix<float> newA = this->A.template toValueType<float>();
+                
+                TopologicalNondeterministicLinearEquationSolver<float> newSolver(newA, this->precision, this->maximalNumberOfIterations, this->relative);
 
-				storm::storage::SparseMatrix<float> new_A = A.template toValueType<float>();
 				std::vector<float> new_x = storm::utility::vector::toValueType<float>(x);
 				std::vector<float> const new_b = storm::utility::vector::toValueType<float>(b);
 
-				tvindles.solveEquationSystem(minimize, new_A, new_x, new_b, nullptr, nullptr);
+				newSolver.solveEquationSystem(minimize, new_x, new_b, nullptr, nullptr);
 
 				for (size_t i = 0, size = new_x.size(); i < size; ++i) {
 					x.at(i) = new_x.at(i);
@@ -86,7 +83,7 @@ namespace storm {
 			}
 
 			// Now, we need to determine the SCCs of the MDP and perform a topological sort.
-			std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = A.getRowGroupIndices();
+			std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = this->A.getRowGroupIndices();
 			
 			// Check if the decomposition is necessary
 #ifdef STORM_HAVE_CUDA
@@ -104,10 +101,7 @@ namespace storm {
 				//std::cout << "Computing the SCC Decomposition took 0ms" << std::endl;
 
 #ifdef STORM_HAVE_CUDA
-				if (!resetCudaDevice()) {
-					LOG4CPLUS_ERROR(logger, "Could not reset CUDA Device, can not use CUDA Equation Solver.");
-					throw storm::exceptions::InvalidStateException() << "Could not reset CUDA Device, can not use CUDA Equation Solver.";
-				}
+                STORM_LOG_THROW(resetCudaDevice(), storm::exceptions::InvalidStateException, "Could not reset CUDA Device, can not use CUDA Equation Solver.");
 
 				std::chrono::high_resolution_clock::time_point calcStartTime = std::chrono::high_resolution_clock::now();
 				bool result = false;
@@ -145,19 +139,19 @@ namespace storm {
 #endif
 			} else {
 				std::chrono::high_resolution_clock::time_point sccStartTime = std::chrono::high_resolution_clock::now();
-				storm::storage::BitVector fullSystem(A.getRowGroupCount(), true);
-				storm::storage::StronglyConnectedComponentDecomposition<ValueType> sccDecomposition(A, fullSystem, false, false);
+				storm::storage::BitVector fullSystem(this->A.getRowGroupCount(), true);
+				storm::storage::StronglyConnectedComponentDecomposition<ValueType> sccDecomposition(this->A, fullSystem, false, false);
 
 				if (sccDecomposition.size() == 0) {
 					LOG4CPLUS_ERROR(logger, "Can not solve given Equation System as the SCC Decomposition returned no SCCs.");
 					throw storm::exceptions::IllegalArgumentException() << "Can not solve given Equation System as the SCC Decomposition returned no SCCs.";
 				}
 
-                storm::storage::SparseMatrix<ValueType> stronglyConnectedComponentsDependencyGraph = sccDecomposition.extractPartitionDependencyGraph(A);
+                storm::storage::SparseMatrix<ValueType> stronglyConnectedComponentsDependencyGraph = sccDecomposition.extractPartitionDependencyGraph(this->A);
 				std::vector<uint_fast64_t> topologicalSort = storm::utility::graph::getTopologicalSort(stronglyConnectedComponentsDependencyGraph);
 
 				// Calculate the optimal distribution of sccs
-				std::vector<std::pair<bool, storm::storage::StateBlock>> optimalSccs = this->getOptimalGroupingFromTopologicalSccDecomposition(sccDecomposition, topologicalSort, A);
+				std::vector<std::pair<bool, storm::storage::StateBlock>> optimalSccs = this->getOptimalGroupingFromTopologicalSccDecomposition(sccDecomposition, topologicalSort, this->A);
 				LOG4CPLUS_INFO(logger, "Optimized SCC Decomposition, originally " << topologicalSort.size() << " SCCs, optimized to " << optimalSccs.size() << " SCCs.");
 
 				std::chrono::high_resolution_clock::time_point sccEndTime = std::chrono::high_resolution_clock::now();
@@ -181,8 +175,8 @@ namespace storm {
 					storm::storage::StateBlock const& scc = sccIndexIt->second;
 
 					// Generate a sub matrix
-					storm::storage::BitVector subMatrixIndices(A.getColumnCount(), scc.cbegin(), scc.cend());
-					storm::storage::SparseMatrix<ValueType> sccSubmatrix = A.getSubmatrix(true, subMatrixIndices, subMatrixIndices);
+					storm::storage::BitVector subMatrixIndices(this->A.getColumnCount(), scc.cbegin(), scc.cend());
+					storm::storage::SparseMatrix<ValueType> sccSubmatrix = this->A.getSubmatrix(true, subMatrixIndices, subMatrixIndices);
 					std::vector<ValueType> sccSubB(sccSubmatrix.getRowCount());
 					storm::utility::vector::selectVectorValues<ValueType>(sccSubB, subMatrixIndices, nondeterministicChoiceIndices, b);
 					std::vector<ValueType> sccSubX(sccSubmatrix.getColumnCount());
@@ -206,7 +200,7 @@ namespace storm {
 						sccSubNondeterministicChoiceIndices.at(outerIndex + 1) = sccSubNondeterministicChoiceIndices.at(outerIndex) + (nondeterministicChoiceIndices[state + 1] - nondeterministicChoiceIndices[state]);
 
 						for (auto rowGroupIt = nondeterministicChoiceIndices[state]; rowGroupIt != nondeterministicChoiceIndices[state + 1]; ++rowGroupIt) {
-							typename storm::storage::SparseMatrix<ValueType>::const_rows row = A.getRow(rowGroupIt);
+							typename storm::storage::SparseMatrix<ValueType>::const_rows row = this->A.getRow(rowGroupIt);
 							for (auto rowIt = row.begin(); rowIt != row.end(); ++rowIt) {
 								if (!subMatrixIndices.get(rowIt->getColumn())) {
 									// This is an outgoing transition of a state in the SCC to a state not included in the SCC
@@ -222,10 +216,7 @@ namespace storm {
 					// For the current SCC, we need to perform value iteration until convergence.
 					if (useGpu) {
 #ifdef STORM_HAVE_CUDA
-						if (!resetCudaDevice()) {
-							LOG4CPLUS_ERROR(logger, "Could not reset CUDA Device, can not use CUDA Equation Solver.");
-							throw storm::exceptions::InvalidStateException() << "Could not reset CUDA Device, can not use CUDA Equation Solver.";
-						}
+                        STORM_LOG_THROW(resetCudaDevice(), storm::exceptions::InvalidStateException, "Could not reset CUDA Device, can not use CUDA-based equation solver.");
 
 						//LOG4CPLUS_INFO(logger, "Device has " << getTotalCudaMemory() << " Bytes of Memory with " << getFreeCudaMemory() << "Bytes free (" << (static_cast<double>(getFreeCudaMemory()) / static_cast<double>(getTotalCudaMemory())) * 100 << "%).");
 						//LOG4CPLUS_INFO(logger, "We will allocate " << (sizeof(uint_fast64_t)* sccSubmatrix.rowIndications.size() + sizeof(uint_fast64_t)* sccSubmatrix.columnsAndValues.size() * 2 + sizeof(double)* sccSubX.size() + sizeof(double)* sccSubX.size() + sizeof(double)* sccSubB.size() + sizeof(double)* sccSubB.size() + sizeof(uint_fast64_t)* sccSubNondeterministicChoiceIndices.size()) << " Bytes.");
@@ -334,7 +325,7 @@ namespace storm {
 
 		template<typename ValueType>
 		std::vector<std::pair<bool, storm::storage::StateBlock>>
-			TopologicalValueIterationNondeterministicLinearEquationSolver<ValueType>::getOptimalGroupingFromTopologicalSccDecomposition(storm::storage::StronglyConnectedComponentDecomposition<ValueType> const& sccDecomposition, std::vector<uint_fast64_t> const& topologicalSort, storm::storage::SparseMatrix<ValueType> const& matrix) const {
+			TopologicalNondeterministicLinearEquationSolver<ValueType>::getOptimalGroupingFromTopologicalSccDecomposition(storm::storage::StronglyConnectedComponentDecomposition<ValueType> const& sccDecomposition, std::vector<uint_fast64_t> const& topologicalSort, storm::storage::SparseMatrix<ValueType> const& matrix) const {
 				std::vector<std::pair<bool, storm::storage::StateBlock>> result;
 #ifdef STORM_HAVE_CUDA
 				// 95% to have a bit of padding
@@ -430,11 +421,11 @@ namespace storm {
 						std::vector<uint_fast64_t> tempGroups;
 						tempGroups.reserve(neededReserveSize);
 						
-						// Copy the first group to make inplace_merge possible
+						// Copy the first group to make inplace_merge possible.
 						storm::storage::StateBlock const& scc_first = sccDecomposition[topologicalSort[startIndex]];
 						tempGroups.insert(tempGroups.cend(), scc_first.cbegin(), scc_first.cend());
 
-						// For set counts <= 80, Inplace Merge is faster
+						// For set counts <= 80, in-place merge is faster.
 						if (((startIndex + 1) + 80) >= topologicalSortSize) {
 							size_t lastSize = 0;
 							for (size_t j = startIndex + 1; j < topologicalSort.size(); ++j) {
@@ -471,7 +462,7 @@ namespace storm {
 		}
 
         // Explicitly instantiate the solver.
-		template class TopologicalValueIterationNondeterministicLinearEquationSolver<double>;
-		template class TopologicalValueIterationNondeterministicLinearEquationSolver<float>;
+		template class TopologicalNondeterministicLinearEquationSolver<double>;
+		template class TopologicalNondeterministicLinearEquationSolver<float>;
     } // namespace solver
 } // namespace storm
diff --git a/src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.h b/src/solver/TopologicalNondeterministicLinearEquationSolver.h
similarity index 88%
rename from src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.h
rename to src/solver/TopologicalNondeterministicLinearEquationSolver.h
index 4ee64b9f8..c8eed7a99 100644
--- a/src/solver/TopologicalValueIterationNondeterministicLinearEquationSolver.h
+++ b/src/solver/TopologicalNondeterministicLinearEquationSolver.h
@@ -19,30 +19,31 @@ namespace storm {
     namespace solver {
         
         /*!
-         * A class that uses SCC Decompositions to solve a linear equation system
+         * A class that uses SCC Decompositions to solve a linear equation system.
          */
         template<class ValueType>
-		class TopologicalValueIterationNondeterministicLinearEquationSolver : public NativeNondeterministicLinearEquationSolver<ValueType> {
+		class TopologicalNondeterministicLinearEquationSolver : public NativeNondeterministicLinearEquationSolver<ValueType> {
         public:
             /*!
              * Constructs a nondeterministic linear equation solver with parameters being set according to the settings
              * object.
+             *
+             * @param A The matrix defining the coefficients of the linear equation system.
              */
-			TopologicalValueIterationNondeterministicLinearEquationSolver();
+			TopologicalNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
             
             /*!
              * Constructs a nondeterminstic linear equation solver with the given parameters.
              *
+             * @param A The matrix defining the coefficients of the linear equation system.
              * @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.
              */
-			TopologicalValueIterationNondeterministicLinearEquationSolver(double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
-            
-            virtual NondeterministicLinearEquationSolver<ValueType>* clone() const override;
+			TopologicalNondeterministicLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
             
-            virtual void solveEquationSystem(bool minimize, storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
+            virtual void solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
 		private:
 
 			bool enableCuda;
diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp
index e18cab3f2..126864c31 100644
--- a/src/utility/solver.cpp
+++ b/src/utility/solver.cpp
@@ -7,6 +7,7 @@
 
 #include "src/solver/NativeNondeterministicLinearEquationSolver.h"
 #include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
+#include "src/solver/TopologicalNondeterministicLinearEquationSolver.h"
 
 #include "src/solver/GurobiLpSolver.h"
 #include "src/solver/GlpkLpSolver.h"
@@ -14,35 +15,72 @@
 namespace storm {
     namespace utility {
         namespace solver {
-            std::unique_ptr<storm::solver::LpSolver> getLpSolver(std::string const& name) {
-                storm::settings::modules::GeneralSettings::LpSolver lpSolver = storm::settings::generalSettings().getLpSolver();
-                switch (lpSolver) {
-                    case storm::settings::modules::GeneralSettings::LpSolver::Gurobi: return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GurobiLpSolver(name));
-                    case storm::settings::modules::GeneralSettings::LpSolver::glpk: return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GlpkLpSolver(name));
+            template<typename ValueType>
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+                storm::settings::modules::GeneralSettings::EquationSolver equationSolver = storm::settings::generalSettings().getEquationSolver();
+                switch (equationSolver) {
+                    case storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
+                    case storm::settings::modules::GeneralSettings::EquationSolver::Native: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix));
                 }
             }
             
             template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> getLinearEquationSolver() {
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+                return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
+            }
+            
+            template<typename ValueType>
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+                return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix));
+            }
+            
+            template<typename ValueType>
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> NondeterministicLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
                 storm::settings::modules::GeneralSettings::EquationSolver equationSolver = storm::settings::generalSettings().getEquationSolver();
                 switch (equationSolver) {
-                    case storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>());
-                    case storm::settings::modules::GeneralSettings::EquationSolver::Native: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>());
+                    case storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx: return std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<ValueType>(matrix));
+                    case storm::settings::modules::GeneralSettings::EquationSolver::Native: return std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>>(new storm::solver::NativeNondeterministicLinearEquationSolver<ValueType>(matrix));
                 }
             }
+
+            template<typename ValueType>
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> GmmxxNondeterministicLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+                return std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<ValueType>(matrix));
+            }
+
+            template<typename ValueType>
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> NativeNondeterministicLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+                return std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>>(new storm::solver::NativeNondeterministicLinearEquationSolver<ValueType>(matrix));
+            }
             
             template<typename ValueType>
-            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> getNondeterministicLinearEquationSolver() {
-                storm::settings::modules::GeneralSettings::EquationSolver equationSolver = storm::settings::generalSettings().getEquationSolver();
-                switch (equationSolver) {
-                    case storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx: return std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<ValueType>());
-                    case storm::settings::modules::GeneralSettings::EquationSolver::Native: return std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>>(new storm::solver::NativeNondeterministicLinearEquationSolver<ValueType>());
+            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> TopologicalNondeterministicLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
+                return std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>>(new storm::solver::TopologicalNondeterministicLinearEquationSolver<ValueType>(matrix));
+            }
+
+            std::unique_ptr<storm::solver::LpSolver> LpSolverFactory::create(std::string const& name) const {
+                storm::settings::modules::GeneralSettings::LpSolver lpSolver = storm::settings::generalSettings().getLpSolver();
+                switch (lpSolver) {
+                    case storm::settings::modules::GeneralSettings::LpSolver::Gurobi: return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GurobiLpSolver(name));
+                    case storm::settings::modules::GeneralSettings::LpSolver::glpk: return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GlpkLpSolver(name));
                 }
             }
             
-            template std::unique_ptr<storm::solver::LinearEquationSolver<double>> getLinearEquationSolver();
+            std::unique_ptr<storm::solver::LpSolver> GlpkLpSolverFactory::create(std::string const& name) const {
+                return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GlpkLpSolver(name));
+            }
 
-            template std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<double>> getNondeterministicLinearEquationSolver();
+            std::unique_ptr<storm::solver::LpSolver> GurobiLpSolverFactory::create(std::string const& name) const {
+                return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GurobiLpSolver(name));
+            }
+            
+            template class LinearEquationSolverFactory<double>;
+            template class GmmxxLinearEquationSolverFactory<double>;
+            template class NativeLinearEquationSolverFactory<double>;
+            template class NondeterministicLinearEquationSolverFactory<double>;
+            template class GmmxxNondeterministicLinearEquationSolverFactory<double>;
+            template class NativeNondeterministicLinearEquationSolverFactory<double>;
+            template class TopologicalNondeterministicLinearEquationSolverFactory<double>;
         }
     }
 }
\ No newline at end of file
diff --git a/src/utility/solver.h b/src/utility/solver.h
index 357a42f0a..137f4e210 100644
--- a/src/utility/solver.h
+++ b/src/utility/solver.h
@@ -10,15 +10,77 @@
 namespace storm {
     namespace utility {
         namespace solver {
+            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;
+            };
             
-            std::unique_ptr<storm::solver::LpSolver> getLpSolver(std::string const& name);
+            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;
+            };
             
             template<typename ValueType>
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> getLinearEquationSolver();
-
+            class GmmxxLinearEquationSolverFactory : public LinearEquationSolverFactory<ValueType> {
+            public:
+                virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            };
+            
+            template<typename ValueType>
+            class NondeterministicLinearEquationSolverFactory {
+            public:
+                /*!
+                 * Creates a new nondeterministic linear equation solver instance with the given matrix.
+                 */
+                virtual std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const;
+            };
+            
             template<typename ValueType>
-            std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> getNondeterministicLinearEquationSolver();
+            class NativeNondeterministicLinearEquationSolverFactory : public NondeterministicLinearEquationSolverFactory<ValueType> {
+            public:
+                virtual std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            };
+            
+            template<typename ValueType>
+            class GmmxxNondeterministicLinearEquationSolverFactory : public NondeterministicLinearEquationSolverFactory<ValueType> {
+            public:
+                virtual std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            };
+            
+            template<typename ValueType>
+            class TopologicalNondeterministicLinearEquationSolverFactory : public NondeterministicLinearEquationSolverFactory<ValueType> {
+            public:
+                virtual std::unique_ptr<storm::solver::NondeterministicLinearEquationSolver<ValueType>> create(storm::storage::SparseMatrix<ValueType> const& matrix) const override;
+            };
 
+            class LpSolverFactory {
+            public:
+                /*!
+                 * Creates a new linear equation solver instance with the given name.
+                 *
+                 * @param name The name of the LP solver.
+                 * @return A pointer to the newly created solver.
+                 */
+                virtual std::unique_ptr<storm::solver::LpSolver> create(std::string const& name) const;
+            };
+            
+            class GlpkLpSolverFactory : public LpSolverFactory {
+            public:
+                virtual std::unique_ptr<storm::solver::LpSolver> create(std::string const& name) const override;
+            };
+            
+            class GurobiLpSolverFactory : public LpSolverFactory {
+            public:
+                virtual std::unique_ptr<storm::solver::LpSolver> create(std::string const& name) const override;
+            };
         }
     }
 }
diff --git a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
index f36b31eca..0ad32ceaf 100644
--- a/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxCtmcCslModelCheckerTest.cpp
@@ -28,7 +28,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Cluster) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::GmmxxLinearEquationSolver<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
     
     // Start checking properties.
     formula = formulaParser.parseFromString("P=? [ F<=100 !\"minimum\"]");
@@ -83,7 +83,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Embedded) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::GmmxxLinearEquationSolver<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
     
     // Start checking properties.
     formula = formulaParser.parseFromString("P=? [ F<=10000 \"down\"]");
@@ -135,7 +135,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Polling) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::GmmxxLinearEquationSolver<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
     
     // Start checking properties.
     formula = formulaParser.parseFromString("P=?[ F<=10 \"target\"]");
@@ -166,7 +166,7 @@ TEST(GmmxxCtmcCslModelCheckerTest, Tandem) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::GmmxxLinearEquationSolver<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
     
     // Start checking properties.
     formula = formulaParser.parseFromString("P=? [ F<=10 \"network_full\" ]");
diff --git a/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
index 0a83e9107..08fb5a29c 100644
--- a/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxDtmcPrctlModelCheckerTest.cpp
@@ -2,7 +2,7 @@
 #include "storm-config.h"
 
 #include "src/logic/Formulas.h"
-#include "src/solver/GmmxxLinearEquationSolver.h"
+#include "src/utility/solver.h"
 #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "src/settings/SettingsManager.h"
@@ -19,7 +19,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, Die) {
 	ASSERT_EQ(dtmc->getNumberOfStates(), 13ull);
 	ASSERT_EQ(dtmc->getNumberOfTransitions(), 20ull);
 
-	storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::GmmxxLinearEquationSolver<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
     
     auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("one");
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
@@ -64,7 +64,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, Crowds) {
 	ASSERT_EQ(8607ull, dtmc->getNumberOfStates());
 	ASSERT_EQ(15113ull, dtmc->getNumberOfTransitions());
 
-    storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::GmmxxLinearEquationSolver<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
 
     auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("observe0Greater1");
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
@@ -100,7 +100,7 @@ TEST(GmmxxDtmcPrctlModelCheckerTest, SynchronousLeader) {
 	ASSERT_EQ(12400ull, dtmc->getNumberOfStates());
 	ASSERT_EQ(16495ull, dtmc->getNumberOfTransitions());
 
-    storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::GmmxxLinearEquationSolver<double>()));
+    storm::modelchecker::SparseDtmcPrctlModelChecker<double> checker(*dtmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxLinearEquationSolverFactory<double>()));
 
     auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
diff --git a/test/functional/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp
index b4b4886a7..0c45cb83d 100644
--- a/test/functional/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/GmmxxMdpPrctlModelCheckerTest.cpp
@@ -2,7 +2,7 @@
 #include "storm-config.h"
 
 #include "src/logic/Formulas.h"
-#include "src/solver/GmmxxNondeterministicLinearEquationSolver.h"
+#include "src/utility/solver.h"
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "src/settings/SettingsManager.h"
@@ -18,7 +18,7 @@ TEST(GmmxxMdpPrctlModelCheckerTest, Dice) {
     ASSERT_EQ(mdp->getNumberOfStates(), 169ull);
     ASSERT_EQ(mdp->getNumberOfTransitions(), 436ull);
     
-    storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::shared_ptr<storm::solver::GmmxxNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
+    storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxNondeterministicLinearEquationSolverFactory<double>()));
     
     auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two");
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
@@ -90,7 +90,7 @@ TEST(GmmxxMdpPrctlModelCheckerTest, Dice) {
     
     std::shared_ptr<storm::models::sparse::Mdp<double>> stateRewardMdp = abstractModel->as<storm::models::sparse::Mdp<double>>();
     
-    storm::modelchecker::SparseMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp, std::shared_ptr<storm::solver::GmmxxNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
+    storm::modelchecker::SparseMdpPrctlModelChecker<double> stateRewardModelChecker(*mdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxNondeterministicLinearEquationSolverFactory<double>()));
     
     labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("done");
     reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
@@ -114,7 +114,7 @@ TEST(GmmxxMdpPrctlModelCheckerTest, Dice) {
     
     std::shared_ptr<storm::models::sparse::Mdp<double>> stateAndTransitionRewardMdp = abstractModel->as<storm::models::sparse::Mdp<double>>();
     
-    storm::modelchecker::SparseMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp, std::shared_ptr<storm::solver::GmmxxNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
+    storm::modelchecker::SparseMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxNondeterministicLinearEquationSolverFactory<double>()));
     
     labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("done");
     reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
@@ -143,7 +143,7 @@ TEST(GmmxxMdpPrctlModelCheckerTest, AsynchronousLeader) {
     ASSERT_EQ(3172ull, mdp->getNumberOfStates());
     ASSERT_EQ(7144ull, mdp->getNumberOfTransitions());
     
-    storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::shared_ptr<storm::solver::GmmxxNondeterministicLinearEquationSolver<double>>(new storm::solver::GmmxxNondeterministicLinearEquationSolver<double>()));
+    storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::GmmxxNondeterministicLinearEquationSolverFactory<double>()));
     
     auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
diff --git a/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp b/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp
index 169510c2a..0c7e4cca7 100644
--- a/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseCtmcCslModelCheckerTest.cpp
@@ -6,7 +6,7 @@
 #include "src/logic/Formulas.h"
 #include "src/builder/ExplicitPrismModelBuilder.h"
 
-#include "src/solver/NativeLinearEquationSolver.h"
+#include "src/utility/solver.h"
 #include "src/modelchecker/csl/SparseCtmcCslModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 
@@ -28,7 +28,7 @@ TEST(SparseCtmcCslModelCheckerTest, Cluster) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::NativeLinearEquationSolver<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
 
     // Start checking properties.
     formula = formulaParser.parseFromString("P=? [ F<=100 !\"minimum\"]");
@@ -83,7 +83,7 @@ TEST(SparseCtmcCslModelCheckerTest, Embedded) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::NativeLinearEquationSolver<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
 
     // Start checking properties.
     formula = formulaParser.parseFromString("P=? [ F<=10000 \"down\"]");
@@ -135,7 +135,7 @@ TEST(SparseCtmcCslModelCheckerTest, Polling) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc);
+    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
     
     // Start checking properties.
     formula = formulaParser.parseFromString("P=?[ F<=10 \"target\"]");
@@ -166,7 +166,7 @@ TEST(SparseCtmcCslModelCheckerTest, Tandem) {
     uint_fast64_t initialState = *ctmc->getInitialStates().begin();
     
     // Create model checker.
-    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::solver::LinearEquationSolver<double>>(new storm::solver::NativeLinearEquationSolver<double>()));
+    storm::modelchecker::SparseCtmcCslModelChecker<double> modelchecker(*ctmc, std::unique_ptr<storm::utility::solver::LinearEquationSolverFactory<double>>(new storm::utility::solver::NativeLinearEquationSolverFactory<double>()));
     
     // Start checking properties.
     formula = formulaParser.parseFromString("P=? [ F<=10 \"network_full\" ]");
diff --git a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
index 3db01e977..2587ceafd 100644
--- a/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/SparseMdpPrctlModelCheckerTest.cpp
@@ -2,7 +2,7 @@
 #include "storm-config.h"
 
 #include "src/logic/Formulas.h"
-#include "src/solver/NativeNondeterministicLinearEquationSolver.h"
+#include "src/utility/solver.h"
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "src/settings/SettingsManager.h"
@@ -18,7 +18,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 169ull);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 436ull);
     
-    storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
+    storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeNondeterministicLinearEquationSolverFactory<double>()));
     
     auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two");
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
@@ -90,7 +90,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	std::shared_ptr<storm::models::sparse::Mdp<double>> stateRewardMdp = abstractModel->as<storm::models::sparse::Mdp<double>>();
     
-    storm::modelchecker::SparseMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
+    storm::modelchecker::SparseMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeNondeterministicLinearEquationSolverFactory<double>()));
 
     labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("done");
     reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
@@ -114,7 +114,7 @@ TEST(SparseMdpPrctlModelCheckerTest, Dice) {
     
 	std::shared_ptr<storm::models::sparse::Mdp<double>> stateAndTransitionRewardMdp = abstractModel->as<storm::models::sparse::Mdp<double>>();
     
-	storm::modelchecker::SparseMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
+	storm::modelchecker::SparseMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeNondeterministicLinearEquationSolverFactory<double>()));
     
     labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("done");
     reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(labelFormula);
@@ -143,7 +143,7 @@ TEST(SparseMdpPrctlModelCheckerTest, AsynchronousLeader) {
 	ASSERT_EQ(3172ull, mdp->getNumberOfStates());
 	ASSERT_EQ(7144ull, mdp->getNumberOfTransitions());
 
-    storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::shared_ptr<storm::solver::NativeNondeterministicLinearEquationSolver<double>>(new storm::solver::NativeNondeterministicLinearEquationSolver<double>()));
+    storm::modelchecker::SparseMdpPrctlModelChecker<double> checker(*mdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::NativeNondeterministicLinearEquationSolverFactory<double>()));
 
     auto labelFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
     auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(labelFormula);
diff --git a/test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp b/test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp
index 471986ce6..932516c34 100644
--- a/test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp
+++ b/test/functional/modelchecker/TopologicalValueIterationMdpPrctlModelCheckerTest.cpp
@@ -3,8 +3,7 @@
 
 
 #include "src/logic/Formulas.h"
-#include "src/solver/NativeNondeterministicLinearEquationSolver.h"
-#include "src/modelchecker/prctl/TopologicalValueIterationMdpPrctlModelChecker.h"
+#include "src/utility/solver.h"
 #include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "src/settings/SettingsManager.h"
@@ -20,7 +19,7 @@ TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, Dice) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 169ull);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 436ull);
     
-	storm::modelchecker::TopologicalValueIterationMdpPrctlModelChecker<double> mc(*mdp);
+    storm::modelchecker::SparseMdpPrctlModelChecker<double> mc(*mdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::TopologicalNondeterministicLinearEquationSolverFactory<double>()));
     
 	//storm::property::prctl::Ap<double>* apFormula = new storm::property::prctl::Ap<double>("two");
 	auto apFormula = std::make_shared<storm::logic::AtomicLabelFormula>("two");
@@ -138,8 +137,7 @@ TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, Dice) {
 	// ------------- state rewards --------------
 	std::shared_ptr<storm::models::sparse::Mdp<double>> stateRewardMdp = storm::parser::AutoParser::parseModel(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.state.rew", "")->as<storm::models::sparse::Mdp<double>>();
     
-	storm::modelchecker::TopologicalValueIterationMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp);
-
+	storm::modelchecker::SparseMdpPrctlModelChecker<double> stateRewardModelChecker(*stateRewardMdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::TopologicalNondeterministicLinearEquationSolverFactory<double>()));
 
 	apFormula = std::make_shared<storm::logic::AtomicLabelFormula>("done");
 	reachabilityRewardFormula = std::make_shared<storm::logic::ReachabilityRewardFormula>(apFormula);
@@ -174,7 +172,7 @@ TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, Dice) {
 	// -------------------------------- state and transition reward ------------------------
 	std::shared_ptr<storm::models::sparse::Mdp<double>> stateAndTransitionRewardMdp = storm::parser::AutoParser::parseModel(STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.tra", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.lab", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.state.rew", STORM_CPP_BASE_PATH "/examples/mdp/two_dice/two_dice.flip.trans.rew")->as<storm::models::sparse::Mdp<double>>();
     
-	storm::modelchecker::TopologicalValueIterationMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp);
+	storm::modelchecker::SparseMdpPrctlModelChecker<double> stateAndTransitionRewardModelChecker(*stateAndTransitionRewardMdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::TopologicalNondeterministicLinearEquationSolverFactory<double>()));
     
 
 	apFormula = std::make_shared<storm::logic::AtomicLabelFormula>("done");
@@ -214,7 +212,7 @@ TEST(TopologicalValueIterationMdpPrctlModelCheckerTest, AsynchronousLeader) {
 	ASSERT_EQ(mdp->getNumberOfStates(), 3172ull);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), 7144ull);
 
-	storm::modelchecker::TopologicalValueIterationMdpPrctlModelChecker<double> mc(*mdp);
+	storm::modelchecker::SparseMdpPrctlModelChecker<double> mc(*mdp, std::unique_ptr<storm::utility::solver::NondeterministicLinearEquationSolverFactory<double>>(new storm::utility::solver::TopologicalNondeterministicLinearEquationSolverFactory<double>()));
 
 	auto apFormula = std::make_shared<storm::logic::AtomicLabelFormula>("elected");
 	auto eventuallyFormula = std::make_shared<storm::logic::EventuallyFormula>(apFormula);
diff --git a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
index a044c4f3b..f4db7c409 100644
--- a/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
+++ b/test/functional/solver/GmmxxLinearEquationSolverTest.cpp
@@ -23,10 +23,10 @@ TEST(GmmxxLinearEquationSolver, SolveWithStandardOptions) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver);
+    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(A));
     
-    storm::solver::GmmxxLinearEquationSolver<double> solver;
-    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::gmmxxEquationSolverSettings().getPrecision());
@@ -51,10 +51,10 @@ TEST(GmmxxLinearEquationSolver, gmres) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50));
+    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(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50);
-    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    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::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::gmmxxEquationSolverSettings().getPrecision());
@@ -79,10 +79,10 @@ TEST(GmmxxLinearEquationSolver, qmr) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Qmr, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None));
+    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(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Qmr, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50);
-    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    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));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::gmmxxEquationSolverSettings().getPrecision());
@@ -107,10 +107,10 @@ TEST(GmmxxLinearEquationSolver, bicgstab) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None));
+    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(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Bicgstab, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None);
-    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    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::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::gmmxxEquationSolverSettings().getPrecision());
@@ -135,10 +135,10 @@ TEST(GmmxxLinearEquationSolver, jacobi) {
     std::vector<double> x(3);
     std::vector<double> b = {11, -16, 1};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Jacobi, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None));
+    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(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Jacobi, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None);
-    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Jacobi, 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::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::gmmxxEquationSolverSettings().getPrecision());
@@ -163,10 +163,10 @@ TEST(GmmxxLinearEquationSolver, gmresilu) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Ilu, true, 50));
+    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(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50);
-    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    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::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::gmmxxEquationSolverSettings().getPrecision());
@@ -191,10 +191,10 @@ TEST(GmmxxLinearEquationSolver, gmresdiag) {
     std::vector<double> x(3);
     std::vector<double> b = {16, -4, -7};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxLinearEquationSolver<double> solver(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::Diagonal, true, 50));
+    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(storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, 1e-6, 10000, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, true, 50);
-    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    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::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::gmmxxEquationSolverSettings().getPrecision());
@@ -218,7 +218,7 @@ TEST(GmmxxLinearEquationSolver, MatrixVectorMultplication) {
     std::vector<double> x(5);
     x[4] = 1;
     
-    storm::solver::GmmxxLinearEquationSolver<double> solver;
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(A, x, nullptr, 4));
+    storm::solver::GmmxxLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(x, nullptr, 4));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::gmmxxEquationSolverSettings().getPrecision());
 }
\ No newline at end of file
diff --git a/test/functional/solver/GmmxxNondeterministicLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxNondeterministicLinearEquationSolverTest.cpp
index e2263d38f..0441ebbbf 100644
--- a/test/functional/solver/GmmxxNondeterministicLinearEquationSolverTest.cpp
+++ b/test/functional/solver/GmmxxNondeterministicLinearEquationSolverTest.cpp
@@ -15,11 +15,11 @@ TEST(GmmxxNondeterministicLinearEquationSolver, SolveWithStandardOptions) {
     std::vector<double> x(1);
     std::vector<double> b = {0.099, 0.5};
         
-    storm::solver::GmmxxNondeterministicLinearEquationSolver<double> solver;
-    ASSERT_NO_THROW(solver.solveEquationSystem(true, A, x, b));
+    storm::solver::GmmxxNondeterministicLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.solveEquationSystem(true, x, b));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     
-    ASSERT_NO_THROW(solver.solveEquationSystem(false, A, x, b));
+    ASSERT_NO_THROW(solver.solveEquationSystem(false, x, b));
     ASSERT_LT(std::abs(x[0] - 0.989991), storm::settings::gmmxxEquationSolverSettings().getPrecision());
 }
 
@@ -41,25 +41,25 @@ TEST(GmmxxNondeterministicLinearEquationSolver, MatrixVectorMultiplication) {
     
     std::vector<double> x = {0, 1, 0};
     
-    ASSERT_NO_THROW(storm::solver::GmmxxNondeterministicLinearEquationSolver<double> solver);
+    ASSERT_NO_THROW(storm::solver::GmmxxNondeterministicLinearEquationSolver<double> solver(A));
     
-    storm::solver::GmmxxNondeterministicLinearEquationSolver<double> solver;
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, A, x, nullptr, 1));
+    storm::solver::GmmxxNondeterministicLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, x, nullptr, 1));
     ASSERT_LT(std::abs(x[0] - 0.099), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, A, x, nullptr, 2));
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, x, nullptr, 2));
     ASSERT_LT(std::abs(x[0] - 0.1881), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, A, x, nullptr, 20));
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, x, nullptr, 20));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, A, x, nullptr, 1));
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, x, nullptr, 1));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::gmmxxEquationSolverSettings().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, A, x, nullptr, 20));
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, x, nullptr, 20));
     ASSERT_LT(std::abs(x[0] - 0.9238082658), storm::settings::gmmxxEquationSolverSettings().getPrecision());
 }
\ No newline at end of file
diff --git a/test/functional/solver/NativeLinearEquationSolverTest.cpp b/test/functional/solver/NativeLinearEquationSolverTest.cpp
index d4eccb561..3d30d6b4f 100644
--- a/test/functional/solver/NativeLinearEquationSolverTest.cpp
+++ b/test/functional/solver/NativeLinearEquationSolverTest.cpp
@@ -23,10 +23,10 @@ TEST(NativeLinearEquationSolver, SolveWithStandardOptions) {
     std::vector<double> x(3);
     std::vector<double> b = {11, -16, 1};
     
-    ASSERT_NO_THROW(storm::solver::NativeLinearEquationSolver<double> solver);
+    ASSERT_NO_THROW(storm::solver::NativeLinearEquationSolver<double> solver(A));
     
-    storm::solver::NativeLinearEquationSolver<double> solver;
-    ASSERT_NO_THROW(solver.solveEquationSystem(A, x, b));
+    storm::solver::NativeLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.solveEquationSystem(x, b));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::nativeEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[1] - 3), storm::settings::nativeEquationSolverSettings().getPrecision());
     ASSERT_LT(std::abs(x[2] - (-1)), storm::settings::nativeEquationSolverSettings().getPrecision());
@@ -50,7 +50,7 @@ TEST(NativeLinearEquationSolver, MatrixVectorMultplication) {
     std::vector<double> x(5);
     x[4] = 1;
     
-    storm::solver::NativeLinearEquationSolver<double> solver;
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(A, x, nullptr, 4));
+    storm::solver::NativeLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(x, nullptr, 4));
     ASSERT_LT(std::abs(x[0] - 1), storm::settings::nativeEquationSolverSettings().getPrecision());
 }
\ No newline at end of file
diff --git a/test/functional/solver/NativeNondeterministicLinearEquationSolverTest.cpp b/test/functional/solver/NativeNondeterministicLinearEquationSolverTest.cpp
index 652c88ee1..55cf5fb2f 100644
--- a/test/functional/solver/NativeNondeterministicLinearEquationSolverTest.cpp
+++ b/test/functional/solver/NativeNondeterministicLinearEquationSolverTest.cpp
@@ -15,11 +15,11 @@ TEST(NativeNondeterministicLinearEquationSolver, SolveWithStandardOptions) {
     std::vector<double> x(1);
     std::vector<double> b = {0.099, 0.5};
         
-    storm::solver::NativeNondeterministicLinearEquationSolver<double> solver;
-    ASSERT_NO_THROW(solver.solveEquationSystem(true, A, x, b));
+    storm::solver::NativeNondeterministicLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.solveEquationSystem(true, x, b));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::nativeEquationSolverSettings().getPrecision());
     
-    ASSERT_NO_THROW(solver.solveEquationSystem(false, A, x, b));
+    ASSERT_NO_THROW(solver.solveEquationSystem(false, x, b));
     ASSERT_LT(std::abs(x[0] - 0.989991), storm::settings::nativeEquationSolverSettings().getPrecision());
 }
 
@@ -41,25 +41,25 @@ TEST(NativeNondeterministicLinearEquationSolver, MatrixVectorMultiplication) {
     
     std::vector<double> x = {0, 1, 0};
     
-    ASSERT_NO_THROW(storm::solver::NativeNondeterministicLinearEquationSolver<double> solver);
+    ASSERT_NO_THROW(storm::solver::NativeNondeterministicLinearEquationSolver<double> solver(A));
     
-    storm::solver::NativeNondeterministicLinearEquationSolver<double> solver;
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, A, x, nullptr, 1));
+    storm::solver::NativeNondeterministicLinearEquationSolver<double> solver(A);
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, x, nullptr, 1));
     ASSERT_LT(std::abs(x[0] - 0.099), storm::settings::nativeEquationSolverSettings().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, A, x, nullptr, 2));
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, x, nullptr, 2));
     ASSERT_LT(std::abs(x[0] - 0.1881), storm::settings::nativeEquationSolverSettings().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, A, x, nullptr, 20));
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(true, x, nullptr, 20));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::nativeEquationSolverSettings().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, A, x, nullptr, 1));
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, x, nullptr, 1));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::nativeEquationSolverSettings().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, A, x, nullptr, 20));
+    ASSERT_NO_THROW(solver.performMatrixVectorMultiplication(false, x, nullptr, 20));
     ASSERT_LT(std::abs(x[0] - 0.9238082658), storm::settings::nativeEquationSolverSettings().getPrecision());
 }
\ No newline at end of file