diff --git a/src/settings/modules/GeneralSettings.cpp b/src/settings/modules/GeneralSettings.cpp
index 7a9ae2971..1274deebf 100644
--- a/src/settings/modules/GeneralSettings.cpp
+++ b/src/settings/modules/GeneralSettings.cpp
@@ -6,7 +6,7 @@
 #include "src/settings/OptionBuilder.h"
 #include "src/settings/ArgumentBuilder.h"
 #include "src/settings/Argument.h"
-
+#include "src/solver/SolverSelectionOptions.h"
 
 #include "src/exceptions/InvalidSettingsException.h"
 
@@ -241,12 +241,12 @@ namespace storm {
                 return this->getOption(timeoutOptionName).getArgumentByName("time").getValueAsUnsignedInteger();
             }
             
-            GeneralSettings::EquationSolver GeneralSettings::getEquationSolver() const {
+            storm::solver::EquationSolverType  GeneralSettings::getEquationSolver() const {
                 std::string equationSolverName = this->getOption(eqSolverOptionName).getArgumentByName("name").getValueAsString();
                 if (equationSolverName == "gmm++") {
-                    return GeneralSettings::EquationSolver::Gmmxx;
+                    return storm::solver::EquationSolverType::Gmmxx;
                 } else if (equationSolverName == "native") {
-                    return GeneralSettings::EquationSolver::Native;
+                    return storm::solver::EquationSolverType::Native;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown equation solver '" << equationSolverName << "'.");
             }
@@ -255,12 +255,12 @@ namespace storm {
                 return this->getOption(eqSolverOptionName).getHasOptionBeenSet();
             }
             
-            GeneralSettings::LpSolver GeneralSettings::getLpSolver() const {
+            storm::solver::LpSolverType GeneralSettings::getLpSolver() const {
                 std::string lpSolverName = this->getOption(lpSolverOptionName).getArgumentByName("name").getValueAsString();
                 if (lpSolverName == "gurobi") {
-                    return GeneralSettings::LpSolver::Gurobi;
+                    return storm::solver::LpSolverType::Gurobi;
                 } else if (lpSolverName == "glpk") {
-                    return GeneralSettings::LpSolver::glpk;
+                    return storm::solver::LpSolverType::Glpk;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown LP solver '" << lpSolverName << "'.");
             }
@@ -297,12 +297,12 @@ namespace storm {
                 return this->getOption(prismCompatibilityOptionName).getHasOptionBeenSet();
             }
 
-			GeneralSettings::MinMaxTechnique GeneralSettings::getMinMaxEquationSolvingTechnique() const {
+			storm::solver::MinMaxTechnique GeneralSettings::getMinMaxEquationSolvingTechnique() const {
 				std::string minMaxEquationSolvingTechnique = this->getOption(minMaxEquationSolvingTechniqueOptionName).getArgumentByName("name").getValueAsString();
 				if (minMaxEquationSolvingTechnique == "valueIteration") {
-					return GeneralSettings::MinMaxTechnique::ValueIteration;
+					return storm::solver::MinMaxTechnique::ValueIteration;
 				} else if (minMaxEquationSolvingTechnique == "policyIteration") {
-					return GeneralSettings::MinMaxTechnique::PolicyIteration;
+					return storm::solver::MinMaxTechnique::PolicyIteration;
 				}
 				STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown min/max equation solving technique '" << minMaxEquationSolvingTechnique << "'.");
 			}
diff --git a/src/settings/modules/GeneralSettings.h b/src/settings/modules/GeneralSettings.h
index f2454a28f..d344528af 100644
--- a/src/settings/modules/GeneralSettings.h
+++ b/src/settings/modules/GeneralSettings.h
@@ -5,6 +5,12 @@
 #include "src/settings/modules/ModuleSettings.h"
 
 namespace storm {
+    namespace solver {
+        enum class EquationSolverType;
+        enum class LpSolverType;
+        enum class MinMaxTechnique;
+    }
+    
     namespace settings {
         namespace modules {
             
@@ -16,14 +22,6 @@ namespace storm {
                 // An enumeration of all engines.
                 enum class Engine { Sparse, Hybrid, Dd };
                 
-                // An enumeration of all available LP solvers.
-                enum class LpSolver { Gurobi, glpk };
-                
-                // An enumeration of all available equation solvers.
-                enum class EquationSolver { Gmmxx, Native };
-
-				// An enumeration of all available Min/Max equation solver techniques.
-				enum class MinMaxTechnique { ValueIteration, PolicyIteration };
                 
                 /*!
                  * Creates a new set of general settings that is managed by the given manager.
@@ -265,7 +263,7 @@ namespace storm {
                  *
                  * @return The selected convergence criterion.
                  */
-                EquationSolver getEquationSolver() const;
+                storm::solver::EquationSolverType getEquationSolver() const;
                 
                 /*!
                  * Retrieves whether a equation solver has been set.
@@ -279,7 +277,7 @@ namespace storm {
                  *
                  * @return The selected LP solver.
                  */
-                LpSolver getLpSolver() const;
+                storm::solver::LpSolverType getLpSolver() const;
                 
                 /*!
                  * Retrieves whether the export-to-dot option was set.
@@ -350,7 +348,7 @@ namespace storm {
 				*
 				* @return The selected min/max equation solving technique.
 				*/
-				MinMaxTechnique getMinMaxEquationSolvingTechnique() const;
+				storm::solver::MinMaxTechnique getMinMaxEquationSolvingTechnique() const;
 
                 
                 bool check() const override;
diff --git a/src/settings/modules/GlpkSettings.cpp b/src/settings/modules/GlpkSettings.cpp
index 7474feb00..a459ddbb8 100644
--- a/src/settings/modules/GlpkSettings.cpp
+++ b/src/settings/modules/GlpkSettings.cpp
@@ -6,6 +6,7 @@
 
 #include "src/settings/SettingsManager.h"
 #include "src/settings/modules/GeneralSettings.h"
+#include "src/solver/SolverSelectionOptions.h"
 
 namespace storm {
     namespace settings {
@@ -34,7 +35,7 @@ namespace storm {
             
             bool GlpkSettings::check() const {
                 if (isOutputSet() || isIntegerToleranceSet()) {
-                    STORM_LOG_WARN_COND(storm::settings::generalSettings().getLpSolver() == storm::settings::modules::GeneralSettings::LpSolver::glpk, "glpk is not selected as the used LP solver, so setting options for glpk has no effect.");
+                    STORM_LOG_WARN_COND(storm::settings::generalSettings().getLpSolver() == storm::solver::LpSolverType::Glpk, "glpk is not selected as the preferred LP solver, so setting options for glpk might have no effect.");
                 }
                 
                 return true;
diff --git a/src/settings/modules/GmmxxEquationSolverSettings.cpp b/src/settings/modules/GmmxxEquationSolverSettings.cpp
index dbaf8832a..00526dbf3 100644
--- a/src/settings/modules/GmmxxEquationSolverSettings.cpp
+++ b/src/settings/modules/GmmxxEquationSolverSettings.cpp
@@ -7,6 +7,7 @@
 
 #include "src/settings/SettingsManager.h"
 #include "src/settings/modules/GeneralSettings.h"
+#include "src/solver/SolverSelectionOptions.h"
 
 namespace storm {
     namespace settings {
@@ -107,7 +108,7 @@ namespace storm {
             bool GmmxxEquationSolverSettings::check() const {
                 bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
                 
-                STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx || !optionsSet, "gmm++ is not selected as the equation solver, so setting options for gmm++ has no effect.");
+                STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::solver::EquationSolverType::Gmmxx || !optionsSet, "gmm++ is not selected as the preferred equation solver, so setting options for gmm++ might have no effect.");
                 
                 return true;
             }
diff --git a/src/settings/modules/GurobiSettings.cpp b/src/settings/modules/GurobiSettings.cpp
index b669a45d3..837c631de 100644
--- a/src/settings/modules/GurobiSettings.cpp
+++ b/src/settings/modules/GurobiSettings.cpp
@@ -6,7 +6,7 @@
 #include "src/settings/Argument.h"
 #include "src/settings/SettingsManager.h"
 #include "src/settings/modules/GeneralSettings.h"
-
+#include "src/solver/SolverSelectionOptions.h"
 namespace storm {
     namespace settings {
         namespace modules {
@@ -46,7 +46,7 @@ namespace storm {
             
             bool GurobiSettings::check() const {
                 if (isOutputSet() || isIntegerToleranceSet() || isNumberOfThreadsSet()) {
-                    STORM_LOG_WARN_COND(storm::settings::generalSettings().getLpSolver() == storm::settings::modules::GeneralSettings::LpSolver::Gurobi, "Gurobi is not selected as the used LP solver, so setting options for Gurobi has no effect.");
+                    STORM_LOG_WARN_COND(storm::settings::generalSettings().getLpSolver() == storm::solver::LpSolverType::Gurobi, "Gurobi is not selected as the preferred LP solver, so setting options for Gurobi might have no effect.");
                 }
                 
                 return true;
diff --git a/src/settings/modules/NativeEquationSolverSettings.cpp b/src/settings/modules/NativeEquationSolverSettings.cpp
index 5584308c5..fb61ddf16 100644
--- a/src/settings/modules/NativeEquationSolverSettings.cpp
+++ b/src/settings/modules/NativeEquationSolverSettings.cpp
@@ -6,7 +6,7 @@
 #include "src/settings/OptionBuilder.h"
 #include "src/settings/ArgumentBuilder.h"
 #include "src/settings/Argument.h"
-
+#include "src/solver/SolverSelectionOptions.h"
 
 namespace storm {
     namespace settings {
@@ -80,7 +80,7 @@ namespace storm {
             bool NativeEquationSolverSettings::check() const {
                 bool optionSet = isLinearEquationSystemTechniqueSet() || isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
                 
-                STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Native || !optionSet, "Native is not selected as the equation solver, so setting options for native has no effect.");
+                STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::solver::EquationSolverType::Native || !optionSet, "Native is not selected as the preferred equation solver, so setting options for native might have no effect.");
                 
                 return true;
             }
diff --git a/src/settings/modules/TopologicalValueIterationEquationSolverSettings.cpp b/src/settings/modules/TopologicalValueIterationEquationSolverSettings.cpp
index 1be47979a..a1a00e283 100644
--- a/src/settings/modules/TopologicalValueIterationEquationSolverSettings.cpp
+++ b/src/settings/modules/TopologicalValueIterationEquationSolverSettings.cpp
@@ -7,6 +7,7 @@
 
 #include "src/settings/SettingsManager.h"
 #include "src/settings/modules/GeneralSettings.h"
+#include "src/solver/SolverSelectionOptions.h"
 
 namespace storm {
     namespace settings {
@@ -54,7 +55,7 @@ namespace storm {
 			bool TopologicalValueIterationEquationSolverSettings::check() const {
                 bool optionsSet = isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet();
                 
-                STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx || !optionsSet, "gmm++ is not selected as the equation solver, so setting options for gmm++ has no effect.");
+                //STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx || !optionsSet, "gmm++ is not selected as the equation solver, so setting options for gmm++ has no effect.");
                 
                 return true;
             }
diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
index a92130b7e..b4357b6d4 100644
--- a/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxMinMaxLinearEquationSolver.cpp
@@ -14,27 +14,24 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        GmmxxMinMaxLinearEquationSolver<ValueType>::GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), stormMatrix(A), rowGroupIndices(A.getRowGroupIndices()) {
-            // Get the settings object to customize solving.
-            storm::settings::modules::GmmxxEquationSolverSettings const& settings = storm::settings::gmmxxEquationSolverSettings();
-			storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings();
-            
-            // Get appropriate settings.
-			precision = settings.getPrecision();            
-            relative = settings.getConvergenceCriterion() == storm::settings::modules::GmmxxEquationSolverSettings::ConvergenceCriterion::Relative;
-			maximalNumberOfIterations = settings.getMaximalIterationCount();
-			useValueIteration = (generalSettings.getMinMaxEquationSolvingTechnique() == storm::settings::modules::GeneralSettings::MinMaxTechnique::ValueIteration);
+        GmmxxMinMaxLinearEquationSolver<ValueType>::GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, MinMaxTechniqueSelection preferredTechnique, bool trackPolicy) : 
+        MinMaxLinearEquationSolver<ValueType>(A, storm::settings::gmmxxEquationSolverSettings().getPrecision(), \
+                                                storm::settings::gmmxxEquationSolverSettings().getConvergenceCriterion() == storm::settings::modules::GmmxxEquationSolverSettings::ConvergenceCriterion::Relative,\
+                                                storm::settings::gmmxxEquationSolverSettings().getMaximalIterationCount(), trackPolicy, preferredTechnique),
+            gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), rowGroupIndices(A.getRowGroupIndices()) {
+           
+			
         }
         
         template<typename ValueType>
-		GmmxxMinMaxLinearEquationSolver<ValueType>::GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool useValueIteration, bool relative) : gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), stormMatrix(A), rowGroupIndices(A.getRowGroupIndices()), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations), useValueIteration(useValueIteration) {
+		GmmxxMinMaxLinearEquationSolver<ValueType>::GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, MinMaxTechniqueSelection tech, bool relative, bool trackPolicy) : MinMaxLinearEquationSolver<ValueType>(A, precision, relative, maximalNumberOfIterations, trackPolicy, tech), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)),  rowGroupIndices(A.getRowGroupIndices())  {
             // Intentionally left empty.
         }
 
         
         template<typename ValueType>
         void GmmxxMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
-			if (useValueIteration) {
+			if (this->useValueIteration) {
 				// 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) {
@@ -56,18 +53,14 @@ namespace storm {
 
 				// Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number
 				// of iterations.
-				while (!converged && iterations < maximalNumberOfIterations) {
+				while (!converged && iterations < this->maximalNumberOfIterations) {
 					// Compute x' = A*x + b.
 					gmm::mult(*gmmxxMatrix, *currentX, *multiplyResult);
 					gmm::add(b, *multiplyResult);
 
 					// Reduce the vector x by applying min/max over all nondeterministic choices.
-					if (minimize) {
-						storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, rowGroupIndices);
-					} else {
-						storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, rowGroupIndices);
-					}
-
+					storm::utility::vector::reduceVectorMinOrMax(minimize, *multiplyResult, *newX, rowGroupIndices);
+					
 					// Determine whether the method converged.
 					converged = storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->precision, this->relative);
 
@@ -99,7 +92,7 @@ namespace storm {
 			} else {
 				// We will use Policy Iteration to solve the given system.
 				// We first guess an initial choice resolution which will be refined after each iteration.
-				std::vector<typename storm::storage::SparseMatrix<ValueType>::index_type> choiceVector(rowGroupIndices.size() - 1);
+				this->policy = std::vector<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
 
 				// Create our own multiplyResult for solving the deterministic sub-instances.
 				std::vector<ValueType> deterministicMultiplyResult(rowGroupIndices.size() - 1);
@@ -126,13 +119,13 @@ namespace storm {
 
 				// Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number
 				// of iterations.
-				while (!converged && iterations < maximalNumberOfIterations) {
+				while (!converged && iterations < this->maximalNumberOfIterations) {
 					// Take the sub-matrix according to the current choices
-					storm::storage::SparseMatrix<ValueType> submatrix = stormMatrix.selectRowsFromRowGroups(choiceVector, true);
+					storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(this->policy, true);
 					submatrix.convertToEquationSystem();
 					
-					GmmxxLinearEquationSolver<ValueType> gmmLinearEquationSolver(submatrix, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, precision, maximalNumberOfIterations, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, relative, 50);
-					storm::utility::vector::selectVectorValues<ValueType>(subB, choiceVector, rowGroupIndices, b);
+					GmmxxLinearEquationSolver<ValueType> gmmLinearEquationSolver(submatrix, storm::solver::GmmxxLinearEquationSolver<double>::SolutionMethod::Gmres, this->precision, this->maximalNumberOfIterations, storm::solver::GmmxxLinearEquationSolver<double>::Preconditioner::None, this->relative, 50);
+					storm::utility::vector::selectVectorValues<ValueType>(subB, this->policy, rowGroupIndices, b);
 
 					// Copy X since we will overwrite it
 					std::copy(currentX->begin(), currentX->end(), newX->begin());
@@ -146,11 +139,8 @@ namespace storm {
 
 					// Reduce the vector x by applying min/max over all nondeterministic choices.
 					// Here, we capture which choice was taken in each state, thereby refining our initial guess.
-					if (minimize) {
-						storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, rowGroupIndices, &choiceVector);
-					} else {
-						storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, rowGroupIndices, &choiceVector);
-					}
+					storm::utility::vector::reduceVectorMinOrMax(minimize, *multiplyResult, *newX, rowGroupIndices, &(this->policy));
+					
 
 					// Determine whether the method converged.
 					converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, static_cast<ValueType>(this->precision), this->relative);
@@ -199,11 +189,8 @@ namespace storm {
                     gmm::add(*b, *multiplyResult);
                 }
                 
-                if (minimize) {
-                    storm::utility::vector::reduceVectorMin(*multiplyResult, x, rowGroupIndices);
-                } else {
-                    storm::utility::vector::reduceVectorMax(*multiplyResult, x, rowGroupIndices);
-                }
+                storm::utility::vector::reduceVectorMinOrMax(minimize, *multiplyResult, x, rowGroupIndices);
+                
             }
             
             if (!multiplyResultMemoryProvided) {
diff --git a/src/solver/GmmxxMinMaxLinearEquationSolver.h b/src/solver/GmmxxMinMaxLinearEquationSolver.h
index 4d151d706..b34f761cb 100644
--- a/src/solver/GmmxxMinMaxLinearEquationSolver.h
+++ b/src/solver/GmmxxMinMaxLinearEquationSolver.h
@@ -20,7 +20,7 @@ namespace storm {
              *
              * @param A The matrix defining the coefficients of the linear equation system.             
              */
-            GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, MinMaxTechniqueSelection preferredTechnique = MinMaxTechniqueSelection::FROMSETTINGS, bool trackPolicy = false);
             
             /*!
              * Constructs a min/max linear equation solver with the given parameters.
@@ -31,7 +31,7 @@ namespace storm {
              * @param relative If set, the relative error rather than the absolute error is considered for convergence
              * detection.
              */
-            GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool useValueIteration, bool relative = true);
+            GmmxxMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, MinMaxTechniqueSelection tech, bool relative, bool trackPolicy = false);
 
             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;
             
@@ -41,23 +41,9 @@ namespace storm {
             // The (gmm++) matrix associated with this equation solver.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix;
 
-			// A reference to the original sparse matrix.
-			storm::storage::SparseMatrix<ValueType> const& stormMatrix;
-            
             // 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;
-            
-            // Sets whether the relative or absolute error is to be considered for convergence detection.
-            bool relative;
-            
-            // The maximal number of iterations to do before iteration is aborted.
-            uint_fast64_t maximalNumberOfIterations;
-
-			// Whether value iteration or policy iteration is to be used.
-			bool useValueIteration;
         };
     } // namespace solver
 } // namespace storm
diff --git a/src/solver/MinMaxLinearEquationSolver.cpp b/src/solver/MinMaxLinearEquationSolver.cpp
index 2aad22bda..680ee7ee0 100644
--- a/src/solver/MinMaxLinearEquationSolver.cpp
+++ b/src/solver/MinMaxLinearEquationSolver.cpp
@@ -1,6 +1,9 @@
 #include "MinMaxLinearEquationSolver.h"
 #include "src/settings/SettingsManager.h"
 #include "src/settings/modules/GeneralSettings.h"
+
+#include "src/utility/macros.h"
+#include "src/exceptions/NotImplementedException.h"
 #include <cstdint>
 
 namespace storm {
@@ -15,5 +18,10 @@ namespace storm {
                 useValueIteration = (prefTech == MinMaxTechniqueSelection::ValueIteration);
             }
         }
+        
+        std::vector<storm::storage::sparse::state_type> AbstractMinMaxLinearEquationSolver::getPolicy() const {
+            STORM_LOG_THROW(!useValueIteration, storm::exceptions::NotImplementedException, "Getting policies after value iteration is not yet supported!");
+            return policy;
+        }
     }
 }
diff --git a/src/solver/MinMaxLinearEquationSolver.h b/src/solver/MinMaxLinearEquationSolver.h
index e9679a6b4..6e6887cd8 100644
--- a/src/solver/MinMaxLinearEquationSolver.h
+++ b/src/solver/MinMaxLinearEquationSolver.h
@@ -2,7 +2,9 @@
 #define STORM_SOLVER_MINMAXLINEAREQUATIONSOLVER_H_
 
 #include <vector>
-
+#include <cstdint>
+#include "SolverSelectionOptions.h"
+#include "src/storage/sparse/StateType.h"
 
 namespace storm {
     namespace storage {
@@ -10,6 +12,38 @@ namespace storm {
     }
     
     namespace solver {
+        /**
+         * Abstract base class which provides value-type independent helpers.
+         */
+        class AbstractMinMaxLinearEquationSolver {
+        
+        public:
+            void setPolicyTracking(bool setToTrue=true);
+            
+            std::vector<storm::storage::sparse::state_type> getPolicy() const;
+            
+        protected:
+            AbstractMinMaxLinearEquationSolver(double precision, bool relativeError, uint_fast64_t maximalIterations, bool trackPolicy, MinMaxTechniqueSelection prefTech);
+             
+            
+            /// The required precision for the iterative methods.
+            double precision;
+            
+            /// Sets whether the relative or absolute error is to be considered for convergence detection.
+            bool relative;
+            
+            /// The maximal number of iterations to do before iteration is aborted.
+            uint_fast64_t maximalNumberOfIterations;
+
+            /// Whether value iteration or policy iteration is to be used.
+            bool useValueIteration;
+            
+            /// Whether we track the policy we generate.
+            bool trackPolicy;
+            
+            /// 
+            mutable std::vector<storm::storage::sparse::state_type> policy;
+        };
         
         /*!
          * A interface that represents an abstract nondeterministic linear equation solver. In addition to solving
@@ -17,7 +51,14 @@ namespace storm {
          * provided.
          */
         template<class ValueType>
-        class MinMaxLinearEquationSolver {
+        class MinMaxLinearEquationSolver : public AbstractMinMaxLinearEquationSolver {
+        protected:
+            MinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& matrix, double precision, bool relativeError, uint_fast64_t maxNrIterations, bool trackPolicy, MinMaxTechniqueSelection prefTech) :
+                AbstractMinMaxLinearEquationSolver(precision, relativeError, maxNrIterations, trackPolicy, prefTech),
+                A(matrix) {
+                // Intentionally left empty.
+            }
+        
         public:
             
             virtual ~MinMaxLinearEquationSolver() {
@@ -59,6 +100,9 @@ namespace storm {
              * @return The result of the repeated matrix-vector multiplication as the content of the vector x.
              */
             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;
+            
+        protected:
+            storm::storage::SparseMatrix<ValueType> const& A;
         };
         
     } // namespace solver
diff --git a/src/solver/NativeMinMaxLinearEquationSolver.cpp b/src/solver/NativeMinMaxLinearEquationSolver.cpp
index 22eb207ad..87ee3b043 100644
--- a/src/solver/NativeMinMaxLinearEquationSolver.cpp
+++ b/src/solver/NativeMinMaxLinearEquationSolver.cpp
@@ -14,30 +14,29 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        NativeMinMaxLinearEquationSolver<ValueType>::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : A(A) {
-            // Get the settings object to customize solving.
-            storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
-			storm::settings::modules::GeneralSettings const& generalSettings = storm::settings::generalSettings();
-            
-            // Get appropriate settings.
-			precision = settings.getPrecision();
-			relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
-            maximalNumberOfIterations = settings.getMaximalIterationCount();
-			useValueIteration = (generalSettings.getMinMaxEquationSolvingTechnique() == storm::settings::modules::GeneralSettings::MinMaxTechnique::ValueIteration);
+        NativeMinMaxLinearEquationSolver<ValueType>::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, MinMaxTechniqueSelection preferredTechnique, bool trackPolicy) : 
+        MinMaxLinearEquationSolver<ValueType>(A, storm::settings::nativeEquationSolverSettings().getPrecision(), \
+            storm::settings::nativeEquationSolverSettings().getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative, \
+            storm::settings::nativeEquationSolverSettings().getMaximalIterationCount(), trackPolicy, preferredTechnique)
+        {
+            // Intentionally left empty.
         }
         
         template<typename ValueType>
-		NativeMinMaxLinearEquationSolver<ValueType>::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool useValueIteration, bool relative) : A(A), precision(precision), relative(relative), maximalNumberOfIterations(maximalNumberOfIterations), useValueIteration(useValueIteration) {
+		NativeMinMaxLinearEquationSolver<ValueType>::NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, MinMaxTechniqueSelection tech, bool relative, bool trackPolicy) : 
+        MinMaxLinearEquationSolver<ValueType>(A, precision, \
+            relative, \
+            maximalNumberOfIterations, trackPolicy, tech) {
             // Intentionally left empty.
         }
         
         template<typename ValueType>
         void NativeMinMaxLinearEquationSolver<ValueType>::solveEquationSystem(bool minimize, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
-			if (useValueIteration) {
+			if (this->useValueIteration) {
 				// 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>(this->A.getRowCount());
 					multiplyResultMemoryProvided = false;
 				}
 				std::vector<ValueType>* currentX = &x;
@@ -54,21 +53,17 @@ namespace storm {
 
 				// Proceed with the iterations as long as the method did not converge or reach the
 				// user-specified maximum number of iterations.
-				while (!converged && iterations < maximalNumberOfIterations) {
+				while (!converged && iterations < this->maximalNumberOfIterations) {
 					// Compute x' = A*x + b.
-					A.multiplyWithVector(*currentX, *multiplyResult);
+					this->A.multiplyWithVector(*currentX, *multiplyResult);
 					storm::utility::vector::addVectors(*multiplyResult, b, *multiplyResult);
 
 					// Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
 					// element of the min/max operator stack.
-					if (minimize) {
-						storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, A.getRowGroupIndices());
-					} else {
-						storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, A.getRowGroupIndices());
-					}
-
+					storm::utility::vector::reduceVectorMinOrMax(minimize, *multiplyResult, *newX, this->A.getRowGroupIndices());
+					
 					// Determine whether the method converged.
-					converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, static_cast<ValueType>(precision), relative);
+					converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, static_cast<ValueType>(this->precision), this->relative);
 
 					// Update environment variables.
 					std::swap(currentX, newX);
@@ -98,11 +93,11 @@ namespace storm {
 			} else {
 				// We will use Policy Iteration to solve the given system.
 				// We first guess an initial choice resolution which will be refined after each iteration.
-				std::vector<typename storm::storage::SparseMatrix<ValueType>::index_type> choiceVector(A.getRowGroupIndices().size() - 1);
+				this->policy = std::vector<storm::storage::sparse::state_type>(this->A.getRowGroupIndices().size() - 1);
 
 				// Create our own multiplyResult for solving the deterministic sub-instances.
-				std::vector<ValueType> deterministicMultiplyResult(A.getRowGroupIndices().size() - 1);
-				std::vector<ValueType> subB(A.getRowGroupIndices().size() - 1);
+				std::vector<ValueType> deterministicMultiplyResult(this->A.getRowGroupIndices().size() - 1);
+				std::vector<ValueType> subB(this->A.getRowGroupIndices().size() - 1);
 
 				// Check whether intermediate storage was provided and create it otherwise.
 				bool multiplyResultMemoryProvided = true;
@@ -126,13 +121,13 @@ namespace storm {
 
 				// Proceed with the iterations as long as the method did not converge or reach the user-specified maximum number
 				// of iterations.
-				while (!converged && iterations < maximalNumberOfIterations) {
+				while (!converged && iterations < this->maximalNumberOfIterations) {
 					// Take the sub-matrix according to the current choices
-					storm::storage::SparseMatrix<ValueType> submatrix = A.selectRowsFromRowGroups(choiceVector, true);
+					storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(this->policy, true);
 					submatrix.convertToEquationSystem();
 
 					NativeLinearEquationSolver<ValueType> nativeLinearEquationSolver(submatrix);
-					storm::utility::vector::selectVectorValues<ValueType>(subB, choiceVector, A.getRowGroupIndices(), b);
+					storm::utility::vector::selectVectorValues<ValueType>(subB, this->policy, this->A.getRowGroupIndices(), b);
 
 					// Copy X since we will overwrite it
 					std::copy(currentX->begin(), currentX->end(), newX->begin());
@@ -141,16 +136,13 @@ namespace storm {
 					nativeLinearEquationSolver.solveEquationSystem(*newX, subB, &deterministicMultiplyResult);
 
 					// Compute x' = A*x + b. This step is necessary to allow the choosing of the optimal policy for the next iteration.
-					A.multiplyWithVector(*newX, *multiplyResult);
+					this->A.multiplyWithVector(*newX, *multiplyResult);
 					storm::utility::vector::addVectors(*multiplyResult, b, *multiplyResult);
 
 					// Reduce the vector x by applying min/max over all nondeterministic choices.
 					// Here, we capture which choice was taken in each state, thereby refining our initial guess.
-					if (minimize) {
-						storm::utility::vector::reduceVectorMin(*multiplyResult, *newX, A.getRowGroupIndices(), &choiceVector);
-					} else {
-						storm::utility::vector::reduceVectorMax(*multiplyResult, *newX, A.getRowGroupIndices(), &choiceVector);
-					}
+					storm::utility::vector::reduceVectorMinOrMax(minimize, *multiplyResult, *newX, this->A.getRowGroupIndices(), &(this->policy));
+					
 
 					// Determine whether the method converged.
 					converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, static_cast<ValueType>(this->precision), this->relative);
@@ -166,6 +158,7 @@ namespace storm {
 				} else {
 					LOG4CPLUS_WARN(logger, "Iterative solver did not converge after " << iterations << " iterations.");
 				}
+                
 
 				// If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
 				// is currently stored in currentX, but x is the output vector.
@@ -189,13 +182,13 @@ namespace storm {
             // 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>(this->A.getRowCount());
                 multiplyResultMemoryProvided = false;
             }
 
             // Now perform matrix-vector multiplication as long as we meet the bound of the formula.
             for (uint_fast64_t i = 0; i < n; ++i) {
-                A.multiplyWithVector(x, *multiplyResult);
+                this->A.multiplyWithVector(x, *multiplyResult);
                 
                 // Add b if it is non-null.
                 if (b != nullptr) {
@@ -204,11 +197,8 @@ namespace storm {
                 
                 // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
                 // element of the min/max operator stack.
-                if (minimize) {
-                    storm::utility::vector::reduceVectorMin(*multiplyResult, x, A.getRowGroupIndices());
-                } else {
-                    storm::utility::vector::reduceVectorMax(*multiplyResult, x, A.getRowGroupIndices());
-                }
+                storm::utility::vector::reduceVectorMinOrMax(minimize, *multiplyResult, x, this->A.getRowGroupIndices());
+                
             }
             
             if (!multiplyResultMemoryProvided) {
diff --git a/src/solver/NativeMinMaxLinearEquationSolver.h b/src/solver/NativeMinMaxLinearEquationSolver.h
index c3fa66b35..5e297096a 100644
--- a/src/solver/NativeMinMaxLinearEquationSolver.h
+++ b/src/solver/NativeMinMaxLinearEquationSolver.h
@@ -19,7 +19,7 @@ namespace storm {
              *
              * @param A The matrix defining the coefficients of the linear equation system.
              */
-            NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, MinMaxTechniqueSelection preferredTechnique = MinMaxTechniqueSelection::FROMSETTINGS, bool trackPolicy = false);
             
             /*!
              * Constructs a min/max linear equation solver with the given parameters.
@@ -30,27 +30,12 @@ namespace storm {
              * @param relative If set, the relative error rather than the absolute error is considered for convergence
              * detection.
              */
-			NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool useValueIteration, bool relative = true);
+            NativeMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, MinMaxTechniqueSelection tech, bool relative = true, bool trackPolicy = false);
             
             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, 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;
-            
-            // Sets whether the relative or absolute error is to be considered for convergence detection.
-            bool relative;
-            
-            // The maximal number of iterations to do before iteration is aborted.
-            uint_fast64_t maximalNumberOfIterations;
-
-			// Whether value iteration or policy iteration is to be used.
-			bool useValueIteration;
         };
     } // namespace solver
 } // namespace storm
diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
index d3d122e12..d31a02d3e 100644
--- a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
+++ b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
@@ -27,19 +27,12 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        TopologicalMinMaxLinearEquationSolver<ValueType>::TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : NativeMinMaxLinearEquationSolver<ValueType>(A) {
+        TopologicalMinMaxLinearEquationSolver<ValueType>::TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A) : 
+        NativeMinMaxLinearEquationSolver<ValueType>(A, storm::settings::topologicalValueIterationEquationSolverSettings().getPrecision(), \
+                storm::settings::topologicalValueIterationEquationSolverSettings().getMaximalIterationCount(), MinMaxTechniqueSelection::ValueIteration, \
+                storm::settings::topologicalValueIterationEquationSolverSettings().getConvergenceCriterion() == storm::settings::modules::TopologicalValueIterationEquationSolverSettings::ConvergenceCriterion::Relative)
+        {
 			// Get the settings object to customize solving.
-			
-			//storm::settings::Settings* settings = storm::settings::Settings::getInstance();
-			auto settings = storm::settings::topologicalValueIterationEquationSolverSettings();
-			// Get appropriate settings.
-			//this->maximalNumberOfIterations = settings->getOptionByLongName("maxiter").getArgument(0).getValueAsUnsignedInteger();
-			//this->precision = settings->getOptionByLongName("precision").getArgument(0).getValueAsDouble();
-			//this->relative = !settings->isSet("absolute");
-			this->maximalNumberOfIterations = settings.getMaximalIterationCount();
-			this->precision = settings.getPrecision();
-			this->relative = (settings.getConvergenceCriterion() == storm::settings::modules::TopologicalValueIterationEquationSolverSettings::ConvergenceCriterion::Relative);
-
 			auto generalSettings = storm::settings::generalSettings();
 			this->enableCuda = generalSettings.isCudaSet();
 #ifdef STORM_HAVE_CUDA
@@ -48,7 +41,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-		TopologicalMinMaxLinearEquationSolver<ValueType>::TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : NativeMinMaxLinearEquationSolver<ValueType>(A, precision, maximalNumberOfIterations, relative) {
+		TopologicalMinMaxLinearEquationSolver<ValueType>::TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative) : NativeMinMaxLinearEquationSolver<ValueType>(A, precision, maximalNumberOfIterations, MinMaxTechniqueSelection::ValueIteration ,relative) {
             // Intentionally left empty.
         }
         
diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.h b/src/solver/TopologicalMinMaxLinearEquationSolver.h
index 4f748672c..be1b844bd 100644
--- a/src/solver/TopologicalMinMaxLinearEquationSolver.h
+++ b/src/solver/TopologicalMinMaxLinearEquationSolver.h
@@ -30,7 +30,7 @@ namespace storm {
              *
              * @param A The matrix defining the coefficients of the linear equation system.
              */
-			TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
+            TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A);
             
             /*!
              * Constructs a min/max linear equation solver with the given parameters.
@@ -41,7 +41,7 @@ namespace storm {
              * @param relative If set, the relative error rather than the absolute error is considered for convergence
              * detection.
              */
-			TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
+            TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision, uint_fast64_t maximalNumberOfIterations, bool relative = true);
             
             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:
diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp
index 53277ec2d..2c9b2fbe8 100644
--- a/src/storage/SparseMatrix.cpp
+++ b/src/storage/SparseMatrix.cpp
@@ -1065,6 +1065,17 @@ namespace storm {
             return sum;
         }
         
+        template<typename ValueType>
+        bool SparseMatrix<ValueType>::isProbabilistic() const {
+            storm::utility::ConstantsComparator<ValueType> comp;
+            for(index_type row = 0; row < this->rowCount; ++row) {
+                if(!comp.isOne(getRowSum(row))) {
+                    return false;
+                }
+            }
+            return true;
+        }
+        
         template<typename ValueType>
         bool SparseMatrix<ValueType>::isSubmatrixOf(SparseMatrix<ValueType> const& matrix) const {
             // Check for matching sizes.
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 1dd16b3de..71fe4078e 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -33,6 +33,9 @@ namespace storm {
         // Forward declare matrix class.
         template<typename T> class SparseMatrix;
         
+        typedef uint_fast64_t SparseMatrixIndexType;
+        
+        
         template<typename IndexType, typename ValueType>
         class MatrixEntry {
         public:
@@ -128,7 +131,7 @@ namespace storm {
         template<typename ValueType>
         class SparseMatrixBuilder {
         public:
-            typedef uint_fast64_t index_type;
+            typedef SparseMatrixIndexType index_type;
             typedef ValueType value_type;
             
             /*!
@@ -300,7 +303,7 @@ namespace storm {
 			friend class storm::solver::TopologicalValueIterationMinMaxLinearEquationSolver<ValueType>;
             friend class SparseMatrixBuilder<ValueType>;
             
-            typedef uint_fast64_t index_type;
+            typedef SparseMatrixIndexType index_type;
             typedef ValueType value_type;
             typedef typename std::vector<MatrixEntry<index_type, value_type>>::iterator iterator;
             typedef typename std::vector<MatrixEntry<index_type, value_type>>::const_iterator const_iterator;
@@ -738,6 +741,10 @@ namespace storm {
              */
             value_type getRowSum(index_type row) const;
             
+            /*!
+             * Checks for each row whether it sums to one.
+             */
+            bool isProbabilistic() const;            
             /*!
              * Checks if the current matrix is a submatrix of the given matrix, where a matrix A is called a submatrix
              * of B if B has no entries in position where A has none. Additionally, the matrices must be of equal size.
@@ -942,6 +949,7 @@ namespace storm {
             
             // A vector indicating the row groups of the matrix.
             std::vector<index_type> rowGroupIndices;
+            
         };
 
     } // namespace storage
diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp
index b58e74cb5..507a72ba3 100644
--- a/src/utility/solver.cpp
+++ b/src/utility/solver.cpp
@@ -24,6 +24,8 @@
 namespace storm {
     namespace utility {
         namespace solver {
+            
+            
             template<storm::dd::DdType Type, typename ValueType>
             std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>> SymbolicLinearEquationSolverFactory<Type, ValueType>::create(storm::dd::Add<Type> const& A, storm::dd::Bdd<Type> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
                 return std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<Type, ValueType>>(new storm::solver::SymbolicLinearEquationSolver<Type, ValueType>(A, allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs));
@@ -36,10 +38,10 @@ namespace storm {
             
             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();
+                storm::solver::EquationSolverType 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));
+                    case storm::solver::EquationSolverType::Gmmxx: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
+                    case storm::solver::EquationSolverType::Native: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::NativeLinearEquationSolver<ValueType>(matrix));
                     default: return std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>>(new storm::solver::GmmxxLinearEquationSolver<ValueType>(matrix));
                 }
             }
@@ -75,10 +77,10 @@ namespace storm {
             
             template<typename ValueType>
             std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> MinMaxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
-                storm::settings::modules::GeneralSettings::EquationSolver equationSolver = storm::settings::generalSettings().getEquationSolver();
+                storm::solver::EquationSolverType equationSolver = storm::settings::generalSettings().getEquationSolver();
                 switch (equationSolver) {
-                    case storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx: return std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>(new storm::solver::GmmxxMinMaxLinearEquationSolver<ValueType>(matrix));
-                    case storm::settings::modules::GeneralSettings::EquationSolver::Native: return std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>(new storm::solver::NativeMinMaxLinearEquationSolver<ValueType>(matrix));
+                    case storm::solver::EquationSolverType::Gmmxx: return std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>(new storm::solver::GmmxxMinMaxLinearEquationSolver<ValueType>(matrix));
+                    case storm::solver::EquationSolverType::Native: return std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>(new storm::solver::NativeMinMaxLinearEquationSolver<ValueType>(matrix));
                     default: return std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>>(new storm::solver::GmmxxMinMaxLinearEquationSolver<ValueType>(matrix));
                 }
             }
@@ -99,10 +101,10 @@ namespace storm {
             }
 
             std::unique_ptr<storm::solver::LpSolver> LpSolverFactory::create(std::string const& name) const {
-                storm::settings::modules::GeneralSettings::LpSolver lpSolver = storm::settings::generalSettings().getLpSolver();
+                storm::solver::LpSolverType 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));
+                    case storm::solver::LpSolverType::Gurobi: return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GurobiLpSolver(name));
+                    case storm::solver::LpSolverType::Glpk: return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GlpkLpSolver(name));
                     default: return std::unique_ptr<storm::solver::LpSolver>(new storm::solver::GurobiLpSolver(name));
                 }
             }
diff --git a/src/utility/solver.h b/src/utility/solver.h
index 661bd90ba..66abbd957 100644
--- a/src/utility/solver.h
+++ b/src/utility/solver.h
@@ -3,6 +3,7 @@
 
 #include "src/solver/NativeLinearEquationSolver.h"
 #include "src/storage/dd/DdType.h"
+#include "src/utility/ExtendSettingEnumWithSelectionField.h"
 
 namespace storm {
     namespace solver {
@@ -22,6 +23,7 @@ namespace storm {
     
     namespace utility {
         namespace solver {
+            
             template<storm::dd::DdType Type, typename ValueType>
             class SymbolicLinearEquationSolverFactory {
             public:
diff --git a/src/utility/vector.h b/src/utility/vector.h
index 634d2bc5a..723550e51 100644
--- a/src/utility/vector.h
+++ b/src/utility/vector.h
@@ -333,7 +333,7 @@ namespace storm {
                 }
 #endif
             }
-            
+                        
             /*!
              * Reduces the given source vector by selecting the smallest element out of each row group.
              *
@@ -360,6 +360,24 @@ namespace storm {
                 reduceVector<T>(source, target, rowGrouping, std::greater<T>(), choices);
             }
             
+            /*!
+             * Reduces the given source vector by selecting either the smallest or the largest out of each row group.
+             *
+             * @param minimize If true, select the smallest, else select the largest.
+             * @param source The source vector which is to be reduced.
+             * @param target The target vector into which a single element from each row group is written.
+             * @param rowGrouping A vector that specifies the begin and end of each group of elements in the source vector.
+             * @param choices If non-null, this vector is used to store the choices made during the selection.
+             */
+            template<class T>
+            void reduceVectorMinOrMax(bool minimize, std::vector<T> const& source, std::vector<T>& target, std::vector<uint_fast64_t> const& rowGrouping, std::vector<uint_fast64_t>* choices = nullptr) {
+                if(minimize) {
+                    reduceVectorMin(source, target, rowGrouping, choices);
+                } else {
+                    reduceVectorMax(source, target, rowGrouping, choices);    
+                }
+            }
+            
             /*!
              * Compares the given elements and determines whether they are equal modulo the given precision. The provided flag
              * additionaly specifies whether the error is computed in relative or absolute terms.
diff --git a/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp
index 5b4174e51..6c659e7b7 100644
--- a/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp
+++ b/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp
@@ -78,7 +78,7 @@ TEST(GmmxxMinMaxLinearEquationSolver, SolveWithPolicyIteration) {
 	std::vector<double> b = { 0.099, 0.5 };
 
 	storm::settings::modules::GmmxxEquationSolverSettings const& settings = storm::settings::gmmxxEquationSolverSettings();
-	storm::solver::GmmxxMinMaxLinearEquationSolver<double> solver(A, settings.getPrecision(), settings.getMaximalIterationCount(), false, settings.getConvergenceCriterion() == storm::settings::modules::GmmxxEquationSolverSettings::ConvergenceCriterion::Relative);
+	storm::solver::GmmxxMinMaxLinearEquationSolver<double> solver(A, settings.getPrecision(), settings.getMaximalIterationCount(), storm::solver::MinMaxTechniqueSelection::PolicyIteration, settings.getConvergenceCriterion() == storm::settings::modules::GmmxxEquationSolverSettings::ConvergenceCriterion::Relative);
 	ASSERT_NO_THROW(solver.solveEquationSystem(true, x, b));
 	ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::gmmxxEquationSolverSettings().getPrecision());
 
diff --git a/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp b/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp
index bb7c281c9..b83423ffe 100644
--- a/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp
+++ b/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp
@@ -79,7 +79,7 @@ TEST(NativeMinMaxLinearEquationSolver, SolveWithPolicyIteration) {
 	std::vector<double> b = { 0.099, 0.5 };
 
 	storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings();
-	storm::solver::NativeMinMaxLinearEquationSolver<double> solver(A, settings.getPrecision(), settings.getMaximalIterationCount(), false, settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative);
+	storm::solver::NativeMinMaxLinearEquationSolver<double> solver(A, settings.getPrecision(), settings.getMaximalIterationCount(), storm::solver::MinMaxTechniqueSelection::PolicyIteration, settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative);
 
 	ASSERT_NO_THROW(solver.solveEquationSystem(true, x, b));
 	ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::nativeEquationSolverSettings().getPrecision());