diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index d9a58c585..4c49a6325 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -89,7 +89,7 @@ namespace storm {
                 }
 
                 // Check for requirements of the solver.
-                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities);
+                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(storm::solver::EquationSystemType::UntilProbabilities);
                 STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver.");
 
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(aProbabilistic);
@@ -376,7 +376,7 @@ namespace storm {
                 std::vector<ValueType> x(numberOfSspStates);
                 
                 // Check for requirements of the solver.
-                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::StochasticShortestPath);
+                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(storm::solver::EquationSystemType::StochasticShortestPath);
                 STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver.");
 
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(sspMatrix);
@@ -584,7 +584,7 @@ namespace storm {
                 std::vector<ValueType> b = probabilisticChoiceRewards;
                 
                 // Check for requirements of the solver.
-                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards);
+                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(storm::solver::EquationSystemType::ReachabilityRewards);
                 STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver.");
 
                 auto solver = minMaxLinearEquationSolverFactory.create(std::move(aProbabilistic));
diff --git a/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
index e89759f81..5c88f9466 100644
--- a/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
@@ -110,7 +110,7 @@ namespace storm {
                         std::vector<ValueType> x(explicitRepresentation.first.getRowGroupCount(), storm::utility::zero<ValueType>());
 
                         // Check for requirements of the solver.
-                        storm::solver::MinMaxLinearEquationSolverRequirements requirements = linearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities, dir);
+                        storm::solver::MinMaxLinearEquationSolverRequirements requirements = linearEquationSolverFactory.getRequirements(storm::solver::EquationSystemType::UntilProbabilities, dir);
                         boost::optional<std::vector<uint64_t>> initialScheduler;
                         if (!requirements.empty()) {
                             if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) {
@@ -351,7 +351,7 @@ namespace storm {
                     // If there are maybe states, we need to solve an equation system.
                     if (!maybeStates.isZero()) {
                         // Check for requirements of the solver this early so we can adapt the maybe states accordingly.
-                        storm::solver::MinMaxLinearEquationSolverRequirements requirements = linearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards, dir);
+                        storm::solver::MinMaxLinearEquationSolverRequirements requirements = linearEquationSolverFactory.getRequirements(storm::solver::EquationSystemType::ReachabilityRewards, dir);
                         bool requireInitialScheduler = false;
                         if (!requirements.empty()) {
                             if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) {
diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
index 5e70c39e9..ae898a14d 100644
--- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
@@ -86,12 +86,12 @@ namespace storm {
             }
             
             template<typename ValueType>
-            std::vector<uint_fast64_t> computeValidSchedulerHint(storm::solver::MinMaxLinearEquationSolverSystemType const& type, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& filterStates, storm::storage::BitVector const& targetStates) {
+            std::vector<uint_fast64_t> computeValidSchedulerHint(storm::solver::EquationSystemType const& type, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& filterStates, storm::storage::BitVector const& targetStates) {
                 storm::storage::Scheduler<ValueType> validScheduler(maybeStates.size());
 
-                if (type == storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities) {
+                if (type == storm::solver::EquationSystemType::UntilProbabilities) {
                     storm::utility::graph::computeSchedulerProbGreater0E(transitionMatrix, backwardTransitions, filterStates, targetStates, validScheduler, boost::none);
-                } else if (type == storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards) {
+                } else if (type == storm::solver::EquationSystemType::ReachabilityRewards) {
                     storm::utility::graph::computeSchedulerProb1E(maybeStates | targetStates, transitionMatrix, backwardTransitions, filterStates, targetStates, validScheduler);
                 } else {
                     STORM_LOG_ASSERT(false, "Unexpected equation system type.");
@@ -198,7 +198,7 @@ namespace storm {
             }
             
             template<typename ValueType>
-            SparseMdpHintType<ValueType> computeHints(storm::solver::MinMaxLinearEquationSolverSystemType const& type, ModelCheckerHint const& hint, storm::OptimizationDirection const& dir, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& targetStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, boost::optional<storm::storage::BitVector> const& selectedChoices = boost::none) {
+            SparseMdpHintType<ValueType> computeHints(storm::solver::EquationSystemType const& type, ModelCheckerHint const& hint, storm::OptimizationDirection const& dir, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, storm::storage::BitVector const& maybeStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& targetStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory, boost::optional<storm::storage::BitVector> const& selectedChoices = boost::none) {
                 SparseMdpHintType<ValueType> result;
 
                 // Check for requirements of the solver.
@@ -216,9 +216,9 @@ namespace storm {
                 extractHintInformationForMaybeStates(result, transitionMatrix, backwardTransitions, maybeStates, selectedChoices, hint, skipEcWithinMaybeStatesCheck);
 
                 result.lowerResultBound = storm::utility::zero<ValueType>();
-                if (type == storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities) {
+                if (type == storm::solver::EquationSystemType::UntilProbabilities) {
                     result.upperResultBound = storm::utility::one<ValueType>();
-                } else if (type == storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards) {
+                } else if (type == storm::solver::EquationSystemType::ReachabilityRewards) {
                     // Intentionally left empty.
                 } else {
                     STORM_LOG_ASSERT(false, "Unexpected equation system type.");
@@ -352,7 +352,7 @@ namespace storm {
                         std::vector<ValueType> b = transitionMatrix.getConstrainedRowGroupSumVector(maybeStates, statesWithProbability1);
                         
                         // Obtain proper hint information either from the provided hint or from requirements of the solver.
-                        SparseMdpHintType<ValueType> hintInformation = computeHints(storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities, hint, goal.direction(), transitionMatrix, backwardTransitions, maybeStates, phiStates, statesWithProbability1, minMaxLinearEquationSolverFactory);
+                        SparseMdpHintType<ValueType> hintInformation = computeHints(storm::solver::EquationSystemType::UntilProbabilities, hint, goal.direction(), transitionMatrix, backwardTransitions, maybeStates, phiStates, statesWithProbability1, minMaxLinearEquationSolverFactory);
                         
                         // Now compute the results for the maybe states.
                         MaybeStateResult<ValueType> resultForMaybeStates = computeValuesForMaybeStates(goal, submatrix, b, produceScheduler, minMaxLinearEquationSolverFactory, hintInformation);
@@ -565,7 +565,7 @@ namespace storm {
                         }
                         
                         // Obtain proper hint information either from the provided hint or from requirements of the solver.
-                        SparseMdpHintType<ValueType> hintInformation = computeHints(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards, hint, goal.direction(), transitionMatrix, backwardTransitions, maybeStates, ~targetStates, targetStates, minMaxLinearEquationSolverFactory, selectedChoices);
+                        SparseMdpHintType<ValueType> hintInformation = computeHints(storm::solver::EquationSystemType::ReachabilityRewards, hint, goal.direction(), transitionMatrix, backwardTransitions, maybeStates, ~targetStates, targetStates, minMaxLinearEquationSolverFactory, selectedChoices);
                         
                         // Now compute the results for the maybe states.
                         MaybeStateResult<ValueType> resultForMaybeStates = computeValuesForMaybeStates(goal, submatrix, b, produceScheduler, minMaxLinearEquationSolverFactory, hintInformation);
@@ -766,7 +766,7 @@ namespace storm {
                 storm::storage::SparseMatrix<ValueType> sspMatrix = sspMatrixBuilder.build(currentChoice, numberOfSspStates, numberOfSspStates);
                 
                 // Check for requirements of the solver.
-                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::StochasticShortestPath);
+                storm::solver::MinMaxLinearEquationSolverRequirements requirements = minMaxLinearEquationSolverFactory.getRequirements(storm::solver::EquationSystemType::StochasticShortestPath);
                 STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver.");
                 
                 std::vector<ValueType> sspResult(numberOfSspStates);
diff --git a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
index a4c418b3c..395018b6e 100644
--- a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
@@ -23,13 +23,13 @@ namespace storm {
         namespace helper {
             
             template<storm::dd::DdType DdType, typename ValueType>
-            storm::dd::Bdd<DdType> computeValidSchedulerHint(storm::solver::MinMaxLinearEquationSolverSystemType const& type, storm::models::symbolic::NondeterministicModel<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& maybeStates, storm::dd::Bdd<DdType> const& targetStates) {
+            storm::dd::Bdd<DdType> computeValidSchedulerHint(storm::solver::EquationSystemType const& type, storm::models::symbolic::NondeterministicModel<DdType, ValueType> const& model, storm::dd::Add<DdType, ValueType> const& transitionMatrix, storm::dd::Bdd<DdType> const& maybeStates, storm::dd::Bdd<DdType> const& targetStates) {
             
                 storm::dd::Bdd<DdType> result;
                 
-                if (type == storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities) {
+                if (type == storm::solver::EquationSystemType::UntilProbabilities) {
                     result = storm::utility::graph::computeSchedulerProbGreater0E(model, transitionMatrix.notZero(), maybeStates, targetStates);
-                } else if (type == storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards) {
+                } else if (type == storm::solver::EquationSystemType::ReachabilityRewards) {
                     result = storm::utility::graph::computeSchedulerProb1E(model, transitionMatrix.notZero(), maybeStates, targetStates, maybeStates || targetStates);
                 }
                 
@@ -82,12 +82,12 @@ namespace storm {
                         std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs());
 
                         // Check requirements of solver.
-                        storm::solver::MinMaxLinearEquationSolverRequirements requirements = solver->getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities, dir);
+                        storm::solver::MinMaxLinearEquationSolverRequirements requirements = solver->getRequirements(storm::solver::EquationSystemType::UntilProbabilities, dir);
                         boost::optional<storm::dd::Bdd<DdType>> initialScheduler;
                         if (!requirements.empty()) {
                             if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) {
                                 STORM_LOG_DEBUG("Computing valid scheduler, because the solver requires it.");
-                                initialScheduler = computeValidSchedulerHint(storm::solver::MinMaxLinearEquationSolverSystemType::UntilProbabilities, model, transitionMatrix, maybeStates, statesWithProbability01.second);
+                                initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::UntilProbabilities, model, transitionMatrix, maybeStates, statesWithProbability01.second);
                                 requirements.set(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler, false);
                             }
                             STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver.");
@@ -241,12 +241,12 @@ namespace storm {
                         std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getIllegalMask() && maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getNondeterminismVariables(), model.getRowColumnMetaVariablePairs());
                         
                         // Check requirements of solver.
-                        storm::solver::MinMaxLinearEquationSolverRequirements requirements = solver->getRequirements(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards, dir);
+                        storm::solver::MinMaxLinearEquationSolverRequirements requirements = solver->getRequirements(storm::solver::EquationSystemType::ReachabilityRewards, dir);
                         boost::optional<storm::dd::Bdd<DdType>> initialScheduler;
                         if (!requirements.empty()) {
                             if (requirements.requires(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler)) {
                                 STORM_LOG_DEBUG("Computing valid scheduler, because the solver requires it.");
-                                initialScheduler = computeValidSchedulerHint(storm::solver::MinMaxLinearEquationSolverSystemType::ReachabilityRewards, model, transitionMatrix, maybeStates, targetStates);
+                                initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::ReachabilityRewards, model, transitionMatrix, maybeStates, targetStates);
                                 requirements.set(storm::solver::MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler, false);
                             }
                             STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver.");
diff --git a/src/storm/solver/MinMaxLinearEquationSolverSystemType.h b/src/storm/solver/EquationSystemType.h
similarity index 76%
rename from src/storm/solver/MinMaxLinearEquationSolverSystemType.h
rename to src/storm/solver/EquationSystemType.h
index e517b7b8f..c68fd3f41 100644
--- a/src/storm/solver/MinMaxLinearEquationSolverSystemType.h
+++ b/src/storm/solver/EquationSystemType.h
@@ -3,7 +3,7 @@
 namespace storm {
     namespace solver {
         
-        enum class MinMaxLinearEquationSolverSystemType {
+        enum class EquationSystemType {
             UntilProbabilities,
             ReachabilityRewards,
             StochasticShortestPath
diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
index f27bd97bf..853f949ef 100644
--- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
@@ -215,16 +215,16 @@ namespace storm {
         }
         
         template<typename ValueType>
-        MinMaxLinearEquationSolverRequirements IterativeMinMaxLinearEquationSolver<ValueType>::getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
+        MinMaxLinearEquationSolverRequirements IterativeMinMaxLinearEquationSolver<ValueType>::getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
             MinMaxLinearEquationSolverRequirements requirements;
             
-            if (equationSystemType == MinMaxLinearEquationSolverSystemType::UntilProbabilities) {
+            if (equationSystemType == EquationSystemType::UntilProbabilities) {
                 if (this->getSettings().getSolutionMethod() == IterativeMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) {
                     if (!direction || direction.get() == OptimizationDirection::Maximize) {
                         requirements.set(MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler);
                     }
                 }
-            } else if (equationSystemType == MinMaxLinearEquationSolverSystemType::ReachabilityRewards) {
+            } else if (equationSystemType == EquationSystemType::ReachabilityRewards) {
                 if (!direction || direction.get() == OptimizationDirection::Minimize) {
                     requirements.set(MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler);
                 }
@@ -235,16 +235,11 @@ namespace storm {
 
         template<typename ValueType>
         bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
-            if(!this->linEqSolverA) {
+            if (!this->linEqSolverA) {
                 this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A);
                 this->linEqSolverA->setCachingEnabled(true);
             }
             
-            if (!this->auxiliaryRowVector) {
-                this->auxiliaryRowVector = std::make_unique<std::vector<ValueType>>(this->A->getRowCount());
-            }
-            std::vector<ValueType>& multiplyResult = *this->auxiliaryRowVector;
-            
             if (!auxiliaryRowGroupVector) {
                 auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
             }
@@ -258,8 +253,12 @@ namespace storm {
                 // Solve the resulting equation system.
                 auto submatrixSolver = this->linearEquationSolverFactory->create(std::move(submatrix));
                 submatrixSolver->setCachingEnabled(true);
-                if (this->lowerBound) { submatrixSolver->setLowerBound(this->lowerBound.get()); }
-                if (this->upperBound) { submatrixSolver->setUpperBound(this->upperBound.get()); }
+                if (this->lowerBound) {
+                    submatrixSolver->setLowerBound(this->lowerBound.get());
+                }
+                if (this->upperBound) {
+                    submatrixSolver->setUpperBound(this->upperBound.get());
+                }
                 submatrixSolver->solveEquations(x, *auxiliaryRowGroupVector);
             }
             
@@ -272,11 +271,8 @@ namespace storm {
             
             Status status = Status::InProgress;
             while (status == Status::InProgress) {
-                // Compute x' = A*x + b.
-                this->linEqSolverA->multiply(*currentX, &b, multiplyResult);
-                
-                // Reduce the vector x' by applying min/max for all non-deterministic choices.
-                storm::utility::vector::reduceVectorMinOrMax(dir, multiplyResult, *newX, this->A->getRowGroupIndices());
+                // Compute x' = min/max(A*x + b).
+                this->linEqSolverA->multiplyAndReduce(dir, this->A->getRowGroupIndices(), *currentX, &b, *newX);
                 
                 // Determine whether the method converged.
                 if (storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) {
@@ -299,14 +295,7 @@ namespace storm {
             
             // If requested, we store the scheduler for retrieval.
             if (this->isTrackSchedulerSet()) {
-                // Due to a custom termination condition, it may be the case that no iterations are performed. In this
-                // case we need to compute x'= A*x+b once.
-                if (iterations==0) {
-                    this->linEqSolverA->multiply(x, &b, multiplyResult);
-                }
-                this->schedulerChoices = std::vector<uint_fast64_t>(this->A->getRowGroupCount());
-                // Reduce the multiplyResult and keep track of the choices made
-                storm::utility::vector::reduceVectorMinOrMax(dir, multiplyResult, x, this->A->getRowGroupIndices(), &this->schedulerChoices.get());
+                this->linEqSolverA->multiplyAndReduce(dir, this->A->getRowGroupIndices(), x, &b, *currentX, &this->schedulerChoices.get());
             }
 
             if (!this->isCachingEnabled()) {
diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h
index 85aeb3843..2b40d0f9a 100644
--- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h
+++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h
@@ -50,7 +50,7 @@ namespace storm {
             ValueType getPrecision() const;
             bool getRelative() const;
             
-            virtual MinMaxLinearEquationSolverRequirements getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const override;
+            virtual MinMaxLinearEquationSolverRequirements getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const override;
             
         private:
             bool solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
diff --git a/src/storm/solver/LinearEquationSolver.cpp b/src/storm/solver/LinearEquationSolver.cpp
index 487b77f89..433c0f419 100644
--- a/src/storm/solver/LinearEquationSolver.cpp
+++ b/src/storm/solver/LinearEquationSolver.cpp
@@ -7,9 +7,14 @@
 #include "storm/solver/EigenLinearEquationSolver.h"
 #include "storm/solver/EliminationLinearEquationSolver.h"
 
+#include "storm/utility/vector.h"
+
 #include "storm/settings/SettingsManager.h"
 #include "storm/settings/modules/CoreSettings.h"
 
+#include "storm/utility/macros.h"
+#include "storm/exceptions/NotSupportedException.h"
+
 namespace storm {
     namespace solver {
         
@@ -20,8 +25,7 @@ namespace storm {
         
         template<typename ValueType>
         void LinearEquationSolver<ValueType>::repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n) const {
-            
-            if(!cachedRowVector) {
+            if (!cachedRowVector) {
                 cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
             }
             
@@ -34,7 +38,6 @@ namespace storm {
             std::vector<ValueType>* currentX = &x;
             std::vector<ValueType>* nextX = cachedRowVector.get();
             
-            
             // Now perform matrix-vector multiplication as long as we meet the bound.
             for (uint_fast64_t i = 0; i < n; ++i) {
                 this->multiply(*currentX, b, *nextX);
@@ -50,11 +53,39 @@ namespace storm {
             // restore the old caching setting
             setCachingEnabled(cachingWasEnabled);
             
-            if(!isCachingEnabled()) {
+            if (!isCachingEnabled()) {
+                clearCache();
+            }
+        }
+        
+        template<typename ValueType>
+        void LinearEquationSolver<ValueType>::multiplyAndReduce(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+            if (!cachedRowVector) {
+                cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
+            }
+            
+            // We enable caching for this. But remember how the old setting was
+            bool cachingWasEnabled = isCachingEnabled();
+            setCachingEnabled(true);
+
+            this->multiply(x, b, *cachedRowVector);
+            storm::utility::vector::reduceVectorMinOrMax(dir, *cachedRowVector, result, rowGroupIndices, choices);
+            
+            // restore the old caching setting
+            setCachingEnabled(cachingWasEnabled);
+            
+            if (!isCachingEnabled()) {
                 clearCache();
             }
         }
         
+#ifdef STORM_HAVE_CARL
+        template<>
+        void LinearEquationSolver<storm::RationalFunction>::multiplyAndReduce(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const* b, std::vector<storm::RationalFunction>& result, std::vector<uint_fast64_t>* choices ) const {
+            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Reducing rational function vector is not supported.");
+        }
+#endif
+        
         template<typename ValueType>
         void LinearEquationSolver<ValueType>::setCachingEnabled(bool value) const {
             if(cachingEnabled && !value) {
diff --git a/src/storm/solver/LinearEquationSolver.h b/src/storm/solver/LinearEquationSolver.h
index 170d625f7..88b8181c5 100644
--- a/src/storm/solver/LinearEquationSolver.h
+++ b/src/storm/solver/LinearEquationSolver.h
@@ -5,6 +5,7 @@
 #include <memory>
 
 #include "storm/solver/AbstractEquationSolver.h"
+#include "storm/solver/OptimizationDirection.h"
 
 #include "storm/storage/SparseMatrix.h"
 
@@ -55,6 +56,22 @@ namespace storm {
              */
             virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const = 0;
             
+            /*!
+             * Performs on matrix-vector multiplication x' = A*x + b and then minimizes/maximizes over the row groups
+             * so that the resulting vector has the size of number of row groups of A.
+             *
+             * @param dir The direction for the reduction step.
+             * @param rowGroupIndices A vector storing the row groups over which to reduce.
+             * @param x The input vector with which to multiply the matrix. Its length must be equal
+             * to the number of columns of A.
+             * @param b If non-null, this vector is added after the multiplication. If given, its length must be equal
+             * to the number of rows of A.
+             * @param result The target vector into which to write the multiplication result. Its length must be equal
+             * to the number of rows of A.
+             * @param choices If given, the choices made in the reduction process are written to this vector.
+             */
+            virtual void multiplyAndReduce(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr) const;
+            
             /*!
              * 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. Note that the
diff --git a/src/storm/solver/MinMaxLinearEquationSolver.cpp b/src/storm/solver/MinMaxLinearEquationSolver.cpp
index ace2aa3db..5fb4fe38b 100644
--- a/src/storm/solver/MinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/MinMaxLinearEquationSolver.cpp
@@ -143,7 +143,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        MinMaxLinearEquationSolverRequirements MinMaxLinearEquationSolver<ValueType>::getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
+        MinMaxLinearEquationSolverRequirements MinMaxLinearEquationSolver<ValueType>::getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
             return MinMaxLinearEquationSolverRequirements();
         }
         
@@ -215,7 +215,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        MinMaxLinearEquationSolverRequirements MinMaxLinearEquationSolverFactory<ValueType>::getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
+        MinMaxLinearEquationSolverRequirements MinMaxLinearEquationSolverFactory<ValueType>::getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
             // Create dummy solver and ask it for requirements.
             std::unique_ptr<MinMaxLinearEquationSolver<ValueType>> solver = this->create();
             return solver->getRequirements(equationSystemType, direction);
diff --git a/src/storm/solver/MinMaxLinearEquationSolver.h b/src/storm/solver/MinMaxLinearEquationSolver.h
index c77417ccd..2d8ee350d 100644
--- a/src/storm/solver/MinMaxLinearEquationSolver.h
+++ b/src/storm/solver/MinMaxLinearEquationSolver.h
@@ -12,7 +12,7 @@
 #include "storm/storage/sparse/StateType.h"
 #include "storm/storage/Scheduler.h"
 #include "storm/solver/OptimizationDirection.h"
-#include "storm/solver/MinMaxLinearEquationSolverSystemType.h"
+#include "storm/solver/EquationSystemType.h"
 #include "storm/solver/MinMaxLinearEquationSolverRequirements.h"
 
 #include "storm/exceptions/InvalidSettingsException.h"
@@ -170,7 +170,7 @@ namespace storm {
              * Retrieves the requirements of this solver for solving equations with the current settings. The requirements
              * are guaranteed to be ordered according to their appearance in the SolverRequirement type.
              */
-            virtual MinMaxLinearEquationSolverRequirements getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const;
+            virtual MinMaxLinearEquationSolverRequirements getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const;
             
             /*!
              * Notifies the solver that the requirements for solving equations have been checked. If this has not been
@@ -233,7 +233,7 @@ namespace storm {
              * Retrieves the requirements of the solver that would be created when calling create() right now. The
              * requirements are guaranteed to be ordered according to their appearance in the SolverRequirement type.
              */
-            MinMaxLinearEquationSolverRequirements getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const;
+            MinMaxLinearEquationSolverRequirements getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const;
             void setRequirementsChecked(bool value = true);
             bool isRequirementsCheckedSet() const;
 
diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp
index ae4824672..ca7cdfb01 100644
--- a/src/storm/solver/NativeLinearEquationSolver.cpp
+++ b/src/storm/solver/NativeLinearEquationSolver.cpp
@@ -223,6 +223,54 @@ namespace storm {
             }
         }
         
+        template<typename ValueType>
+        void NativeLinearEquationSolver<ValueType>::multiplyAndReduce(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+            
+            // If the vector and the result are aliases, we need and temporary vector.
+            std::vector<ValueType>* target;
+            std::vector<ValueType> temporary;
+            if (&x == &result) {
+                STORM_LOG_WARN("Vectors are aliased. Using temporary, may be slow.");
+                temporary = std::vector<ValueType>(x.size());
+                target = &temporary;
+            } else {
+                target = &result;
+            }
+            
+            for (uint64_t rowGroup = 0; rowGroup < rowGroupIndices.size() - 1; ++rowGroup) {
+                uint64_t row = rowGroupIndices[rowGroup];
+
+                (*target)[rowGroup] = b ? (*b)[row] : storm::utility::zero<ValueType>();
+                for (auto const& entry : this->A->getRow(row)) {
+                    (*target)[rowGroup] += entry.getValue() * x[entry.getColumn()];
+                }
+                ++row;
+                
+                for (; row < rowGroupIndices[rowGroup + 1]; ++row) {
+                    ValueType newValue = b ? (*b)[row] : storm::utility::zero<ValueType>();
+                    for (auto const& entry : this->A->getRow(row)) {
+                        newValue += entry.getValue() * x[entry.getColumn()];
+                    }
+                    
+                    if (dir == OptimizationDirection::Minimize && newValue < result[rowGroup]) {
+                        (*target)[rowGroup] = newValue;
+                        if (choices) {
+                            (*choices)[rowGroup] = row - rowGroupIndices[rowGroup];
+                        }
+                    } else if (dir == OptimizationDirection::Maximize && newValue > result[rowGroup]) {
+                        (*target)[rowGroup] = newValue;
+                        if (choices) {
+                            (*choices)[rowGroup] = row - rowGroupIndices[rowGroup];
+                        }
+                    }
+                }
+            }
+            
+            if (!temporary.empty()) {
+                std::swap(temporary, result);
+            }
+        }
+        
         template<typename ValueType>
         void NativeLinearEquationSolver<ValueType>::setSettings(NativeLinearEquationSolverSettings<ValueType> const& newSettings) {
             settings = newSettings;
diff --git a/src/storm/solver/NativeLinearEquationSolver.h b/src/storm/solver/NativeLinearEquationSolver.h
index 2a978cdbc..b5de30ffc 100644
--- a/src/storm/solver/NativeLinearEquationSolver.h
+++ b/src/storm/solver/NativeLinearEquationSolver.h
@@ -51,7 +51,8 @@ namespace storm {
             
             virtual bool solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
             virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
-
+            virtual void multiplyAndReduce(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr) const override;
+            
             void setSettings(NativeLinearEquationSolverSettings<ValueType> const& newSettings);
             NativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
 
diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
index 18fa0fcd3..b4449c45b 100644
--- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
@@ -257,16 +257,16 @@ namespace storm {
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        MinMaxLinearEquationSolverRequirements SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
+        MinMaxLinearEquationSolverRequirements SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
             MinMaxLinearEquationSolverRequirements requirements;
             
-            if (equationSystemType == MinMaxLinearEquationSolverSystemType::UntilProbabilities) {
+            if (equationSystemType == EquationSystemType::UntilProbabilities) {
                 if (this->getSettings().getSolutionMethod() == SymbolicMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) {
                     if (!direction || direction.get() == OptimizationDirection::Maximize) {
                         requirements.set(MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler);
                     }
                 }
-            } else if (equationSystemType == MinMaxLinearEquationSolverSystemType::ReachabilityRewards) {
+            } else if (equationSystemType == EquationSystemType::ReachabilityRewards) {
                 if (!direction || direction.get() == OptimizationDirection::Minimize) {
                     requirements.set(MinMaxLinearEquationSolverRequirements::Element::ValidInitialScheduler);
                 }
@@ -291,7 +291,7 @@ namespace storm {
         }
 
         template<storm::dd::DdType DdType, typename ValueType>
-        MinMaxLinearEquationSolverRequirements SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType>::getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
+        MinMaxLinearEquationSolverRequirements SymbolicMinMaxLinearEquationSolverFactory<DdType, ValueType>::getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction) const {
             std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> solver = this->create();
             return solver->getRequirements(equationSystemType, direction);
         }
diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h
index b0f70ac39..b99e72cd7 100644
--- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h
+++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h
@@ -10,7 +10,7 @@
 
 #include "storm/solver/SymbolicLinearEquationSolver.h"
 
-#include "storm/solver/MinMaxLinearEquationSolverSystemType.h"
+#include "storm/solver/EquationSystemType.h"
 #include "storm/solver/MinMaxLinearEquationSolverRequirements.h"
 
 #include "storm/storage/expressions/Variable.h"
@@ -125,7 +125,7 @@ namespace storm {
             /*!
              * Retrieves the requirements of the solver.
              */
-            virtual MinMaxLinearEquationSolverRequirements getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const;
+            virtual MinMaxLinearEquationSolverRequirements getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const;
             
             /*!
              * Notifies the solver that the requirements for solving equations have been checked. If this has not been
@@ -189,7 +189,7 @@ namespace storm {
              * Retrieves the requirements of the solver that would be created when calling create() right now. The
              * requirements are guaranteed to be ordered according to their appearance in the SolverRequirement type.
              */
-            MinMaxLinearEquationSolverRequirements getRequirements(MinMaxLinearEquationSolverSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const;
+            MinMaxLinearEquationSolverRequirements getRequirements(EquationSystemType const& equationSystemType, boost::optional<storm::solver::OptimizationDirection> const& direction = boost::none) const;
             
         private:
             virtual std::unique_ptr<storm::solver::SymbolicMinMaxLinearEquationSolver<DdType, ValueType>> create() const = 0;