From f2e581b3dfa91f2d0eb2e627176f1c3817c84fcd Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Wed, 11 Oct 2017 15:59:36 +0200
Subject: [PATCH] rational search for symbolic linear equation solvers

---
 .../3rdparty/sylvan/src/sylvan_mtbdd_storm.h  |   4 +-
 .../src/sylvan_storm_rational_function.h      |   2 +-
 .../sylvan/src/sylvan_storm_rational_number.c |  18 +-
 .../sylvan/src/sylvan_storm_rational_number.h |  14 +-
 .../prctl/helper/SymbolicDtmcPrctlHelper.cpp  |  10 +-
 .../prctl/helper/SymbolicMdpPrctlHelper.cpp   |   4 +
 .../SymbolicQuantitativeCheckResult.cpp       |   6 -
 ...ymbolicEliminationLinearEquationSolver.cpp | 117 +++++----
 .../SymbolicEliminationLinearEquationSolver.h |  15 +-
 src/storm/solver/SymbolicEquationSolver.cpp   |  94 +++++++
 src/storm/solver/SymbolicEquationSolver.h     |  54 ++++
 .../solver/SymbolicLinearEquationSolver.cpp   |  83 ++++--
 .../solver/SymbolicLinearEquationSolver.h     |  48 ++--
 .../SymbolicMinMaxLinearEquationSolver.cpp    |  35 +--
 .../SymbolicMinMaxLinearEquationSolver.h      |  12 +-
 .../SymbolicNativeLinearEquationSolver.cpp    | 238 +++++++++++++++++-
 .../SymbolicNativeLinearEquationSolver.h      |  64 ++++-
 .../dd/bisimulation/QuotientExtractor.cpp     |   6 +-
 src/storm/utility/constants.cpp               |   5 +
 19 files changed, 674 insertions(+), 155 deletions(-)
 create mode 100644 src/storm/solver/SymbolicEquationSolver.cpp
 create mode 100644 src/storm/solver/SymbolicEquationSolver.h

diff --git a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h b/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h
index 449a98c3e..d249c5564 100644
--- a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h
+++ b/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h
@@ -164,8 +164,10 @@ TASK_DECL_3(BDD, mtbdd_min_abstract_representative, MTBDD, MTBDD, uint32_t);
 TASK_DECL_3(BDD, mtbdd_max_abstract_representative, MTBDD, MTBDD, uint32_t);
 #define mtbdd_max_abstract_representative(a, vars) (CALL(mtbdd_max_abstract_representative, a, vars, 0))
 
+// A version of unary apply that performs no caching. This is needed of the argument is actually used by the unary operation,
+// but may not be used to identify cache entries, for example if the argument is a pointer.
 TASK_DECL_3(MTBDD, mtbdd_uapply_nocache, MTBDD, mtbdd_uapply_op, size_t);
-#define mtbdd_uapply_nocache(dd, op, param) CALL(mtbdd_uapply_nocache, dd, op, param)
+#define mtbdd_uapply_nocache(dd, op, param) (CALL(mtbdd_uapply_nocache, dd, op, param))
 
 #ifdef __cplusplus
 }
diff --git a/resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h b/resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h
index 888a80d51..b6fdfeed2 100644
--- a/resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h
+++ b/resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h
@@ -74,7 +74,7 @@ TASK_DECL_3(MTBDD, sylvan_storm_rational_function_and_exists, MTBDD, MTBDD, MTBD
 #define sylvan_storm_rational_function_abstract_max(dd, v) mtbdd_abstract(dd, v, TASK(sylvan_storm_rational_function_abstract_op_max))
     
 TASK_DECL_2(MTBDD, sylvan_storm_rational_function_op_to_double, MTBDD, size_t)
-#define sylvan_storm_rational_function_to_double(a) mtbdd_uapply_nocache(a, TASK(sylvan_storm_rational_function_op_to_double), 0)
+#define sylvan_storm_rational_function_to_double(a) mtbdd_uapply(a, TASK(sylvan_storm_rational_function_op_to_double), 0)
     
 TASK_DECL_2(MTBDD, sylvan_storm_rational_function_op_threshold, MTBDD, size_t*)
 TASK_DECL_2(MTBDD, sylvan_storm_rational_function_op_strict_threshold, MTBDD, size_t*)
diff --git a/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c b/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c
index 22f9ef053..261156ba7 100644
--- a/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c
+++ b/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c
@@ -583,8 +583,8 @@ TASK_IMPL_3(MTBDD, sylvan_storm_rational_number_and_exists, MTBDD, a, MTBDD, b,
     return result;
 }
 
-TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, a, size_t*, svalue) {
-    storm_rational_number_ptr value = (storm_rational_number_ptr)svalue;
+TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, a, size_t, svalue) {
+    storm_rational_number_ptr value = (storm_rational_number_ptr)(void*)svalue;
     
     if (mtbdd_isleaf(a)) {
         storm_rational_number_ptr ma = mtbdd_getstorm_rational_number_ptr(a);
@@ -594,8 +594,8 @@ TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, a, size_t*,
     return mtbdd_invalid;
 }
 
-TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, a, size_t*, svalue) {
-    storm_rational_number_ptr value = (storm_rational_number_ptr)svalue;
+TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, a, size_t, svalue) {
+    storm_rational_number_ptr value = (storm_rational_number_ptr)(void*)svalue;
 
     if (mtbdd_isleaf(a)) {
         storm_rational_number_ptr ma = mtbdd_getstorm_rational_number_ptr(a);
@@ -606,16 +606,6 @@ TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, a, s
     return mtbdd_invalid;
 }
 
-TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_threshold, MTBDD, dd, storm_rational_number_ptr, value)
-{
-    return mtbdd_uapply(dd, TASK(sylvan_storm_rational_number_op_threshold), (size_t*)value);
-}
-
-TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_strict_threshold, MTBDD, dd, storm_rational_number_ptr, value)
-{
-    return mtbdd_uapply(dd, TASK(sylvan_storm_rational_number_op_strict_threshold), (size_t*)value);
-}
-
 TASK_IMPL_1(MTBDD, sylvan_storm_rational_number_minimum, MTBDD, a) {
     /* Check terminal case */
     if (a == mtbdd_false) return mtbdd_false;
diff --git a/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h b/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h
index 227a9831d..070f03e1a 100644
--- a/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h
+++ b/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h
@@ -74,16 +74,12 @@ TASK_DECL_3(MTBDD, sylvan_storm_rational_number_and_exists, MTBDD, MTBDD, MTBDD)
 #define sylvan_storm_rational_number_abstract_max(dd, v) mtbdd_abstract(dd, v, TASK(sylvan_storm_rational_number_abstract_op_max))
     
 TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_to_double, MTBDD, size_t)
-#define sylvan_storm_rational_number_to_double(a) mtbdd_uapply_nocache(a, TASK(sylvan_storm_rational_number_op_to_double), 0)
+#define sylvan_storm_rational_number_to_double(a) mtbdd_uapply(a, TASK(sylvan_storm_rational_number_op_to_double), 0)
     
-TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, size_t*)
-TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, size_t*)
-
-TASK_DECL_2(MTBDD, sylvan_storm_rational_number_threshold, MTBDD, storm_rational_number_ptr);
-#define sylvan_storm_rational_number_threshold(dd, value) CALL(sylvan_storm_rational_number_strict_threshold, dd, value)
-    
-TASK_DECL_2(MTBDD, sylvan_storm_rational_number_strict_threshold, MTBDD, storm_rational_number_ptr);
-#define sylvan_storm_rational_number_strict_threshold(dd, value) CALL(sylvan_storm_rational_number_strict_threshold, dd, value)
+TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, size_t)
+TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, size_t)
+#define sylvan_storm_rational_number_threshold(dd, value) mtbdd_uapply_nocache(dd, TASK(sylvan_storm_rational_number_op_threshold), (size_t)(void*)value)
+#define sylvan_storm_rational_number_strict_threshold(dd, value) mtbdd_uapply_nocache(dd, TASK(sylvan_storm_rational_number_op_strict_threshold), (size_t)(void*)value)
     
 TASK_DECL_3(MTBDD, sylvan_storm_rational_number_equal_norm_d, MTBDD, MTBDD, storm_rational_number_ptr);
 #define sylvan_storm_rational_number_equal_norm_d(a, b, epsilon) CALL(sylvan_storm_rational_number_equal_norm_d, a, b, epsilon)
diff --git a/src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp
index 86f624459..ffc602aaa 100644
--- a/src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp
@@ -52,10 +52,13 @@ namespace storm {
                         // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
                         // for solving the equation system (i.e. compute (I-A)).
                         submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
-                        submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                        if (linearEquationSolverFactory.getEquationProblemFormat() == storm::solver::LinearEquationSolverProblemFormat::EquationSystem) {
+                            submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                        }
                         
                         // Solve the equation system.
                         std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                        solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
                         storm::dd::Add<DdType, ValueType> result = solver->solveEquations(model.getManager().template getAddZero<ValueType>(), subvector);
                         
                         return statesWithProbability01.second.template toAdd<ValueType>() + result;
@@ -168,10 +171,13 @@ namespace storm {
                         // Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
                         // for solving the equation system (i.e. compute (I-A)).
                         submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
-                        submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                        if (linearEquationSolverFactory.getEquationProblemFormat() == storm::solver::LinearEquationSolverProblemFormat::EquationSystem) {
+                            submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+                        }
                         
                         // Solve the equation system.
                         std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
+                        solver->setLowerBound(storm::utility::zero<ValueType>());
                         storm::dd::Add<DdType, ValueType> result = solver->solveEquations(model.getManager().template getAddZero<ValueType>(), subvector);
                         
                         return infinityStates.ite(model.getManager().getConstant(storm::utility::infinity<ValueType>()), result);
diff --git a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
index bca6bec0f..518f50acc 100644
--- a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
+++ b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
@@ -87,11 +87,13 @@ namespace storm {
                                 initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::UntilProbabilities, model, transitionMatrix, maybeStates, statesWithProbability01.second);
                                 requirements.clearValidInitialScheduler();
                             }
+                            requirements.clearBounds();
                             STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver.");
                         }
                         if (initialScheduler) {
                             solver->setInitialScheduler(initialScheduler.get());
                         }
+                        solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
                         solver->setRequirementsChecked();
 
                         storm::dd::Add<DdType, ValueType> result = solver->solveEquations(dir, model.getManager().template getAddZero<ValueType>(), subvector);
@@ -244,11 +246,13 @@ namespace storm {
                                 initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::ReachabilityRewards, model, transitionMatrix, maybeStates, targetStates);
                                 requirements.clearValidInitialScheduler();
                             }
+                            requirements.clearLowerBounds();
                             STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver.");
                         }
                         if (initialScheduler) {
                             solver->setInitialScheduler(initialScheduler.get());
                         }
+                        solver->setLowerBound(storm::utility::zero<ValueType>());
                         solver->setRequirementsChecked();
 
                         storm::dd::Add<DdType, ValueType> result = solver->solveEquations(dir, model.getManager().template getAddZero<ValueType>(), subvector);
diff --git a/src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp b/src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
index 6f5597562..b73190bab 100644
--- a/src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
+++ b/src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
@@ -28,13 +28,7 @@ namespace storm {
             if (comparisonType == storm::logic::ComparisonType::Less) {
                 states &= values.less(bound);
             } else if (comparisonType == storm::logic::ComparisonType::LessEqual) {
-                std::cout << "states zero: " << (values.equals(values.getDdManager().template getAddZero<ValueType>()) && reachableStates).getNonZeroCount() << std::endl;
-                std::cout << "states one: " << (values.equals(values.getDdManager().template getAddOne<ValueType>()) && reachableStates).getNonZeroCount() << std::endl;
-                std::cout << "states before: " << states.getNonZeroCount() << std::endl;
-                values.exportToDot("vals.dot");
-                std::cout << "total: " <<((values.equals(values.getDdManager().template getAddOne<ValueType>()) && states) || (values.lessOrEqual(bound) && states)).getNonZeroCount() << std::endl;
                 states &= values.lessOrEqual(bound);
-                std::cout << "states after: " << states.getNonZeroCount() << std::endl;
             } else if (comparisonType == storm::logic::ComparisonType::Greater) {
                 states &= values.greater(bound);
             } else if (comparisonType == storm::logic::ComparisonType::GreaterEqual) {
diff --git a/src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp b/src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp
index 1cdc67743..7443ecdcf 100644
--- a/src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp
@@ -10,6 +10,11 @@
 namespace storm {
     namespace solver {
      
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicEliminationLinearEquationSolver<DdType, ValueType>::SymbolicEliminationLinearEquationSolver() : SymbolicLinearEquationSolver<DdType, ValueType>() {
+            // Intentionally left empty.
+        }
+        
         template<storm::dd::DdType DdType, typename ValueType>
         SymbolicEliminationLinearEquationSolver<DdType, ValueType>::SymbolicEliminationLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> 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) : SymbolicEliminationLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
             this->setMatrix(A);
@@ -17,6 +22,64 @@ namespace storm {
         
         template<storm::dd::DdType DdType, typename ValueType>
         SymbolicEliminationLinearEquationSolver<DdType, ValueType>::SymbolicEliminationLinearEquationSolver(storm::dd::Bdd<DdType> 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) : SymbolicLinearEquationSolver<DdType, ValueType>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
+            this->createInternalData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType, ValueType> SymbolicEliminationLinearEquationSolver<DdType, ValueType>::solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            storm::dd::DdManager<DdType>& ddManager = x.getDdManager();
+            
+            // Build diagonal BDD over new meta variables.
+            storm::dd::Bdd<DdType> diagonal = storm::utility::dd::getRowColumnDiagonal(ddManager, this->rowColumnMetaVariablePairs);
+            diagonal &= this->getAllRows();
+            diagonal = diagonal.swapVariables(this->oldNewMetaVariablePairs);
+
+            storm::dd::Add<DdType, ValueType> rowsAdd = this->getAllRows().swapVariables(rowRowMetaVariablePairs).template toAdd<ValueType>();
+            storm::dd::Add<DdType, ValueType> diagonalAdd = diagonal.template toAdd<ValueType>();
+            
+            // Move the matrix to the new meta variables.
+            storm::dd::Add<DdType, ValueType> matrix = this->A.swapVariables(oldNewMetaVariablePairs);
+            
+            // Initialize solution over the new meta variables.
+            storm::dd::Add<DdType, ValueType> solution = b.swapVariables(oldNewMetaVariablePairs);
+            
+            // As long as there are transitions, we eliminate them.
+            uint64_t iterations = 0;
+            while (!matrix.isZero()) {
+                // Determine inverse loop probabilies.
+                storm::dd::Add<DdType, ValueType> inverseLoopProbabilities = rowsAdd / (rowsAdd - (diagonalAdd * matrix).sumAbstract(newColumnVariables));
+                
+                // Scale all transitions with the inverse loop probabilities.
+                matrix *= inverseLoopProbabilities;
+                
+                solution *= inverseLoopProbabilities;
+                
+                // Delete diagonal elements, i.e. remove self-loops.
+                matrix = diagonal.ite(ddManager.template getAddZero<ValueType>(), matrix);
+            
+                // Update the one-step probabilities.
+                solution += (matrix * solution.swapVariables(newRowColumnMetaVariablePairs)).sumAbstract(newColumnVariables);
+                
+                // Shortcut all transitions by eliminating one intermediate step.
+                matrix = matrix.multiplyMatrix(matrix.permuteVariables(shiftMetaVariablePairs), newColumnVariables);
+                matrix = matrix.swapVariables(columnHelperMetaVariablePairs);
+                
+                ++iterations;
+                STORM_LOG_TRACE("Completed iteration " << iterations << " of elimination process.");
+            }
+            
+            STORM_LOG_INFO("Elimination completed in " << iterations << " iterations.");
+            
+            return solution.swapVariables(rowRowMetaVariablePairs);
+        }
+
+        template<storm::dd::DdType DdType, typename ValueType>
+        LinearEquationSolverProblemFormat SymbolicEliminationLinearEquationSolver<DdType, ValueType>::getEquationProblemFormat() const {
+            return LinearEquationSolverProblemFormat::FixedPointSystem;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        void SymbolicEliminationLinearEquationSolver<DdType, ValueType>::createInternalData(storm::dd::Bdd<DdType> 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) {
             storm::dd::DdManager<DdType>& ddManager = allRows.getDdManager();
             
             // Create triple-layered meta variables for all original meta variables. We will use them later in the elimination process.
@@ -61,56 +124,18 @@ namespace storm {
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        storm::dd::Add<DdType, ValueType> SymbolicEliminationLinearEquationSolver<DdType, ValueType>::solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
-            storm::dd::DdManager<DdType>& ddManager = x.getDdManager();
-            
-            // Build diagonal BDD over new meta variables.
-            storm::dd::Bdd<DdType> diagonal = storm::utility::dd::getRowColumnDiagonal(ddManager, this->rowColumnMetaVariablePairs);
-            diagonal &= this->allRows;
-            diagonal = diagonal.swapVariables(this->oldNewMetaVariablePairs);
-
-            storm::dd::Add<DdType, ValueType> rowsAdd = this->allRows.swapVariables(rowRowMetaVariablePairs).template toAdd<ValueType>();
-            storm::dd::Add<DdType, ValueType> diagonalAdd = diagonal.template toAdd<ValueType>();
-            
-            // Revert the conversion to an equation system and move it to the new meta variables.
-            storm::dd::Add<DdType, ValueType> matrix = diagonalAdd - this->A.swapVariables(oldNewMetaVariablePairs);
+        void SymbolicEliminationLinearEquationSolver<DdType, ValueType>::setData(storm::dd::Bdd<DdType> 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) {
             
-            // Initialize solution over the new meta variables.
-            storm::dd::Add<DdType, ValueType> solution = b.swapVariables(oldNewMetaVariablePairs);
+            // Call superclass function.
+            SymbolicLinearEquationSolver<DdType, ValueType>::setData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
             
-            // As long as there are transitions, we eliminate them.
-            uint64_t iterations = 0;
-            while (!matrix.isZero()) {
-                // Determine inverse loop probabilies.
-                storm::dd::Add<DdType, ValueType> inverseLoopProbabilities = rowsAdd / (rowsAdd - (diagonalAdd * matrix).sumAbstract(newColumnVariables));
-                
-                // Scale all transitions with the inverse loop probabilities.
-                matrix *= inverseLoopProbabilities;
-                
-                solution *= inverseLoopProbabilities;
-                
-                // Delete diagonal elements, i.e. remove self-loops.
-                matrix = diagonal.ite(ddManager.template getAddZero<ValueType>(), matrix);
-            
-                // Update the one-step probabilities.
-                solution += (matrix * solution.swapVariables(newRowColumnMetaVariablePairs)).sumAbstract(newColumnVariables);
-                
-                // Shortcut all transitions by eliminating one intermediate step.
-                matrix = matrix.multiplyMatrix(matrix.permuteVariables(shiftMetaVariablePairs), newColumnVariables);
-                matrix = matrix.swapVariables(columnHelperMetaVariablePairs);
-                
-                ++iterations;
-                STORM_LOG_TRACE("Completed iteration " << iterations << " of elimination process.");
-            }
-            
-            STORM_LOG_INFO("Elimination completed in " << iterations << " iterations.");
-            
-            return solution.swapVariables(rowRowMetaVariablePairs);
+            // Now create new variables as needed.
+            this->createInternalData(this->getAllRows(), this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs);
         }
-
+        
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicEliminationLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Bdd<DdType> 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::make_unique<SymbolicEliminationLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicEliminationLinearEquationSolverFactory<DdType, ValueType>::create() const {
+            return std::make_unique<SymbolicEliminationLinearEquationSolver<DdType, ValueType>>();
         }
             
         template<storm::dd::DdType DdType, typename ValueType>
diff --git a/src/storm/solver/SymbolicEliminationLinearEquationSolver.h b/src/storm/solver/SymbolicEliminationLinearEquationSolver.h
index 45b0f0b12..2e1ce5a4b 100644
--- a/src/storm/solver/SymbolicEliminationLinearEquationSolver.h
+++ b/src/storm/solver/SymbolicEliminationLinearEquationSolver.h
@@ -14,13 +14,26 @@ namespace storm {
         template<storm::dd::DdType DdType, typename ValueType = double>
         class SymbolicEliminationLinearEquationSolver : public SymbolicLinearEquationSolver<DdType, ValueType> {
         public:
+            /*!
+             * Constructs a symbolic linear equation solver.
+             *
+             * @param settings The settings to use.
+             */
+            SymbolicEliminationLinearEquationSolver();
+            
             SymbolicEliminationLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> 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);
 
             SymbolicEliminationLinearEquationSolver(storm::dd::Bdd<DdType> 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);
 
             virtual storm::dd::Add<DdType, ValueType> solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const override;
             
+            virtual void setData(storm::dd::Bdd<DdType> 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) override;
+
+            virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override;
+            
         private:
+            void createInternalData(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
+            
             std::vector<std::vector<storm::expressions::Variable>> oldToNewMapping;
             std::set<storm::expressions::Variable> newRowVariables;
             std::set<storm::expressions::Variable> newColumnVariables;
@@ -38,7 +51,7 @@ namespace storm {
         public:
             using SymbolicLinearEquationSolverFactory<DdType, ValueType>::create;
             
-            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> 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 override;
+            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create() const override;
             
             SymbolicEliminationLinearEquationSolverSettings<ValueType>& getSettings();
             SymbolicEliminationLinearEquationSolverSettings<ValueType> const& getSettings() const;
diff --git a/src/storm/solver/SymbolicEquationSolver.cpp b/src/storm/solver/SymbolicEquationSolver.cpp
new file mode 100644
index 000000000..a76f94873
--- /dev/null
+++ b/src/storm/solver/SymbolicEquationSolver.cpp
@@ -0,0 +1,94 @@
+#include "storm/solver/SymbolicEquationSolver.h"
+
+#include "storm/adapters/RationalNumberAdapter.h"
+#include "storm/adapters/RationalFunctionAdapter.h"
+
+#include "storm/utility/macros.h"
+#include "storm/exceptions/UnmetRequirementException.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicEquationSolver<DdType, ValueType>::SymbolicEquationSolver(storm::dd::Bdd<DdType> const& allRows) : allRows(allRows) {
+            // Intentionally left empty.
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::DdManager<DdType>& SymbolicEquationSolver<DdType, ValueType>::getDdManager() const {
+            return this->allRows.getDdManager();
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        void SymbolicEquationSolver<DdType, ValueType>::setAllRows(storm::dd::Bdd<DdType> const& allRows) {
+            this->allRows = allRows;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Bdd<DdType> const& SymbolicEquationSolver<DdType, ValueType>::getAllRows() const {
+            return this->allRows;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        void SymbolicEquationSolver<DdType, ValueType>::setLowerBounds(storm::dd::Add<DdType, ValueType> const& lowerBounds) {
+            this->lowerBounds = lowerBounds;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        void SymbolicEquationSolver<DdType, ValueType>::setLowerBound(ValueType const& lowerBound) {
+            this->lowerBound = lowerBound;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        void SymbolicEquationSolver<DdType, ValueType>::setUpperBounds(storm::dd::Add<DdType, ValueType> const& upperBounds) {
+            this->upperBounds = upperBounds;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        void SymbolicEquationSolver<DdType, ValueType>::setUpperBound(ValueType const& upperBound) {
+            this->upperBound = upperBound;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        void SymbolicEquationSolver<DdType, ValueType>::setBounds(ValueType const& lowerBound, ValueType const& upperBound) {
+            setLowerBound(lowerBound);
+            setUpperBound(upperBound);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        void SymbolicEquationSolver<DdType, ValueType>::setBounds(storm::dd::Add<DdType, ValueType> const& lowerBounds, storm::dd::Add<DdType, ValueType> const& upperBounds) {
+            setLowerBounds(lowerBounds);
+            setUpperBounds(upperBounds);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType, ValueType> SymbolicEquationSolver<DdType, ValueType>::getLowerBounds() const {
+            STORM_LOG_THROW(lowerBound || lowerBounds, storm::exceptions::UnmetRequirementException, "Requiring lower bounds, but did not get any.");
+            if (lowerBounds) {
+                return lowerBounds.get();
+            } else {
+                return this->allRows.ite(this->allRows.getDdManager().getConstant(lowerBound.get()), this->allRows.getDdManager().template getAddZero<ValueType>());
+            }
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType, ValueType> SymbolicEquationSolver<DdType, ValueType>::getUpperBounds() const {
+            STORM_LOG_THROW(upperBound || upperBounds, storm::exceptions::UnmetRequirementException, "Requiring upper bounds, but did not get any.");
+            if (upperBounds) {
+                return upperBounds.get();
+            } else {
+                return this->allRows.ite(this->allRows.getDdManager().getConstant(upperBound.get()), this->allRows.getDdManager().template getAddZero<ValueType>());
+            }
+        }
+        
+        template class SymbolicEquationSolver<storm::dd::DdType::CUDD, double>;
+        template class SymbolicEquationSolver<storm::dd::DdType::Sylvan, double>;
+        
+#ifdef STORM_HAVE_CARL
+        template class SymbolicEquationSolver<storm::dd::DdType::CUDD, storm::RationalNumber>;
+        template class SymbolicEquationSolver<storm::dd::DdType::Sylvan, storm::RationalNumber>;
+
+        template class SymbolicEquationSolver<storm::dd::DdType::Sylvan, storm::RationalFunction>;
+#endif
+    }
+}
diff --git a/src/storm/solver/SymbolicEquationSolver.h b/src/storm/solver/SymbolicEquationSolver.h
new file mode 100644
index 000000000..2e6322cfc
--- /dev/null
+++ b/src/storm/solver/SymbolicEquationSolver.h
@@ -0,0 +1,54 @@
+#pragma once
+
+#include "storm/storage/dd/DdType.h"
+
+#include "storm/storage/dd/DdManager.h"
+#include "storm/storage/dd/Bdd.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<storm::dd::DdType DdType, typename ValueType = double>
+        class SymbolicEquationSolver {
+        public:
+            SymbolicEquationSolver() = default;
+            SymbolicEquationSolver(storm::dd::Bdd<DdType> const& allRows);
+
+            void setLowerBounds(storm::dd::Add<DdType, ValueType> const& lowerBounds);
+            void setLowerBound(ValueType const& lowerBound);
+            void setUpperBounds(storm::dd::Add<DdType, ValueType> const& upperBounds);
+            void setUpperBound(ValueType const& lowerBound);
+            void setBounds(ValueType const& lowerBound, ValueType const& upperBound);
+            void setBounds(storm::dd::Add<DdType, ValueType> const& lowerBounds, storm::dd::Add<DdType, ValueType> const& upperBounds);
+            
+            /*!
+             * Retrieves a vector of lower bounds for all values (if any lower bounds are known).
+             */
+            storm::dd::Add<DdType, ValueType> getLowerBounds() const;
+            
+            /*!
+             * Retrieves a vector of upper bounds for all values (if any lower bounds are known).
+             */
+            storm::dd::Add<DdType, ValueType> getUpperBounds() const;
+
+        protected:
+            storm::dd::DdManager<DdType>& getDdManager() const;
+
+            void setAllRows(storm::dd::Bdd<DdType> const& allRows);
+            storm::dd::Bdd<DdType> const& getAllRows() const;
+
+            // The relevant rows to this equation solver.
+            storm::dd::Bdd<DdType> allRows;
+
+        private:
+            // Lower bounds (if given).
+            boost::optional<storm::dd::Add<DdType, ValueType>> lowerBounds;
+            boost::optional<ValueType> lowerBound;
+            
+            // Upper bounds (if given).
+            boost::optional<storm::dd::Add<DdType, ValueType>> upperBounds;
+            boost::optional<ValueType> upperBound;
+        };
+        
+    }
+}
diff --git a/src/storm/solver/SymbolicLinearEquationSolver.cpp b/src/storm/solver/SymbolicLinearEquationSolver.cpp
index 41620bbf8..fb6e6059e 100644
--- a/src/storm/solver/SymbolicLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicLinearEquationSolver.cpp
@@ -13,19 +13,25 @@
 #include "storm/settings/modules/CoreSettings.h"
 
 #include "storm/utility/macros.h"
+#include "storm/exceptions/UnmetRequirementException.h"
 
 #include "storm/adapters/RationalFunctionAdapter.h"
 
 namespace storm {
     namespace solver {
 
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver() {
+            // Intentionally left empty.
+        }
+        
         template<storm::dd::DdType DdType, typename ValueType>
         SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> 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) : SymbolicLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
             this->setMatrix(A);
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Bdd<DdType> 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) : allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) {
+        SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Bdd<DdType> 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) : SymbolicEquationSolver<DdType, ValueType>(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) {
             // Intentionally left empty.
         }
         
@@ -45,16 +51,30 @@ namespace storm {
             return xCopy;
         }
         
+        template<storm::dd::DdType DdType, typename ValueType>
+        LinearEquationSolverProblemFormat SymbolicLinearEquationSolver<DdType, ValueType>::getEquationProblemFormat() const {
+            return LinearEquationSolverProblemFormat::EquationSystem;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        LinearEquationSolverRequirements SymbolicLinearEquationSolver<DdType, ValueType>::getRequirements() const {
+            // Return empty requirements by default.
+            return LinearEquationSolverRequirements();
+        }
+        
         template<storm::dd::DdType DdType, typename ValueType>
         void SymbolicLinearEquationSolver<DdType, ValueType>::setMatrix(storm::dd::Add<DdType, ValueType> const& newA) {
             this->A = newA;
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        storm::dd::DdManager<DdType>& SymbolicLinearEquationSolver<DdType, ValueType>::getDdManager() const {
-            return this->allRows.getDdManager();
+        void SymbolicLinearEquationSolver<DdType, ValueType>::setData(storm::dd::Bdd<DdType> 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) {
+            this->setAllRows(allRows);
+            this->rowMetaVariables = rowMetaVariables;
+            this->columnMetaVariables = columnMetaVariables;
+            this->rowColumnMetaVariablePairs = rowColumnMetaVariablePairs;
         }
-        
+                
         template<storm::dd::DdType DdType, typename ValueType>
         std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> 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 {
             std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = this->create(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
@@ -63,55 +83,68 @@ namespace storm {
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> GeneralSymbolicLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Bdd<DdType> 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 {
+        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Bdd<DdType> 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 {
+            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = this->create();
+            solver->setData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+            return solver;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        LinearEquationSolverProblemFormat SymbolicLinearEquationSolverFactory<DdType, ValueType>::getEquationProblemFormat() const {
+            return this->create()->getEquationProblemFormat();
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        LinearEquationSolverRequirements SymbolicLinearEquationSolverFactory<DdType, ValueType>::getRequirements() const {
+            return this->create()->getRequirements();
+        }
+
+        template<storm::dd::DdType DdType, typename ValueType>
+        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> GeneralSymbolicLinearEquationSolverFactory<DdType, ValueType>::create() const {
             storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::CoreSettings>().getEquationSolver();
             switch (equationSolver) {
-                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, ValueType>>();
                     break;
-                case storm::solver::EquationSolverType::Native: return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                case storm::solver::EquationSolverType::Native: return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, ValueType>>();
                     break;
                 default:
-                    STORM_LOG_WARN("The selected equation solver is not available in the dd engine. Falling back to native solver.");
-                    return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                    STORM_LOG_INFO("The selected equation solver is not available in the dd engine. Falling back to native solver.");
+                    return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, ValueType>>();
             }
         }
         
         template<storm::dd::DdType DdType>
-        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalNumber>> GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalNumber>::create(storm::dd::Bdd<DdType> 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 {
+        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalNumber>> GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalNumber>::create() const {
 
             auto const& coreSettings = storm::settings::getModule<storm::settings::modules::CoreSettings>();
             storm::solver::EquationSolverType equationSolver = coreSettings.getEquationSolver();
-            if (coreSettings.isEquationSolverSetFromDefaultValue() && equationSolver != storm::solver::EquationSolverType::Elimination) {
-                STORM_LOG_WARN("Selecting the elimination solver to guarantee exact results. If you want to override this, please explicitly specify a different equation solver.");
-                equationSolver = storm::solver::EquationSolverType::Elimination;
+            if (coreSettings.isEquationSolverSetFromDefaultValue() && equationSolver != storm::solver::EquationSolverType::Native && equationSolver != storm::solver::EquationSolverType::Elimination) {
+                STORM_LOG_INFO("Selecting the native solver to provide a method that guarantees exact results. If you want to override this, please explicitly specify a different equation solver.");
+                equationSolver = storm::solver::EquationSolverType::Native;
             }
             
-            if (equationSolver != storm::solver::EquationSolverType::Elimination) {
-                STORM_LOG_WARN("The chosen equation solver does not guarantee precise results despite using exact arithmetic. Consider using the elimination solver instead.");
-            }
-
             switch (equationSolver) {
-                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalNumber>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalNumber>>();
                     break;
                 case storm::solver::EquationSolverType::Native:
-                    return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, storm::RationalNumber>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                    return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, storm::RationalNumber>>();
                     break;
                 default:
-                    STORM_LOG_WARN("The selected equation solver is not available in the dd engine. Falling back to elimination solver.");
-                    return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalNumber>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                    STORM_LOG_INFO("The selected equation solver is not available in the dd engine. Falling back to elimination solver.");
+                    return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalNumber>>();
             }
         }
         
         template<storm::dd::DdType DdType>
-        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalFunction>> GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalFunction>::create(storm::dd::Bdd<DdType> 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 {
+        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalFunction>> GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalFunction>::create() const {
             
             storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::CoreSettings>().getEquationSolver();
             switch (equationSolver) {
-                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalFunction>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalFunction>>();
                     break;
                 default:
-                    STORM_LOG_WARN("The selected equation solver is not available in the DD setting. Falling back to elimination solver.");
-                    return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalFunction>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
+                    STORM_LOG_INFO("The selected equation solver is not available in the DD setting. Falling back to elimination solver.");
+                    return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalFunction>>();
             }            
         }
         
diff --git a/src/storm/solver/SymbolicLinearEquationSolver.h b/src/storm/solver/SymbolicLinearEquationSolver.h
index 9b3c3b490..d8d2e00ca 100644
--- a/src/storm/solver/SymbolicLinearEquationSolver.h
+++ b/src/storm/solver/SymbolicLinearEquationSolver.h
@@ -8,6 +8,10 @@
 #include "storm/storage/dd/DdManager.h"
 #include "storm/storage/dd/DdType.h"
 
+#include "storm/solver/SymbolicEquationSolver.h"
+#include "storm/solver/LinearEquationSolverProblemFormat.h"
+#include "storm/solver/LinearEquationSolverRequirements.h"
+
 #include "storm/adapters/RationalFunctionAdapter.h"
 
 namespace storm {
@@ -17,24 +21,19 @@ namespace storm {
         
         template<storm::dd::DdType Type>
         class Bdd;
-        
     }
     
     namespace solver {
         
-        template<storm::dd::DdType DdType, typename ValueType = double>
-        class SymbolicLinearEquationSolverSettings {
-        public:
-            // Currently empty.
-        };
-        
         /*!
          * An interface that represents an abstract symbolic linear equation solver. In addition to solving a system of
          * linear equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
          */
         template<storm::dd::DdType DdType, typename ValueType = double>
-        class SymbolicLinearEquationSolver {
+        class SymbolicLinearEquationSolver : public SymbolicEquationSolver<DdType, ValueType> {
         public:
+            SymbolicLinearEquationSolver();
+            
             /*!
              * Constructs a symbolic linear equation solver with the given meta variable sets and pairs.
              *
@@ -85,17 +84,25 @@ namespace storm {
              */
             virtual storm::dd::Add<DdType, ValueType> multiply(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const* b = nullptr, uint_fast64_t n = 1) const;
             
+            /*!
+             * Retrieves the format in which this solver expects to solve equations. If the solver expects the equation
+             * system format, it solves Ax = b. If it it expects a fixed point format, it solves Ax + b = x.
+             */
+            virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const;
+            
+            /*!
+             * Retrieves the requirements of the solver under the current settings. Note that these requirements only
+             * apply to solving linear equations and not to the matrix vector multiplications.
+             */
+            virtual LinearEquationSolverRequirements getRequirements() const;
+            
             void setMatrix(storm::dd::Add<DdType, ValueType> const& newA);
+            virtual void setData(storm::dd::Bdd<DdType> 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);
             
         protected:
-            storm::dd::DdManager<DdType>& getDdManager() const;
-            
             // The matrix defining the coefficients of the linear equation system.
             storm::dd::Add<DdType, ValueType> A;
             
-            // A BDD characterizing all rows of the equation system.
-            storm::dd::Bdd<DdType> allRows;
-            
             // The row variables.
             std::set<storm::expressions::Variable> rowMetaVariables;
             
@@ -103,15 +110,20 @@ namespace storm {
             std::set<storm::expressions::Variable> columnMetaVariables;
             
             // The pairs of meta variables used for renaming.
-            std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs;
+            std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> rowColumnMetaVariablePairs;
         };
         
         template<storm::dd::DdType DdType, typename ValueType>
         class SymbolicLinearEquationSolverFactory {
         public:
-            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> 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 = 0;
+            std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> 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;
             
             std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> 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;
+
+            LinearEquationSolverProblemFormat getEquationProblemFormat() const;
+            LinearEquationSolverRequirements getRequirements() const;
+
+            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create() const = 0;
         };
         
         template<storm::dd::DdType DdType, typename ValueType>
@@ -119,7 +131,7 @@ namespace storm {
         public:
             using SymbolicLinearEquationSolverFactory<DdType, ValueType>::create;
             
-            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> 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;
+            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create() const override;
         };
 
         template<storm::dd::DdType DdType>
@@ -127,7 +139,7 @@ namespace storm {
         public:
             using SymbolicLinearEquationSolverFactory<DdType, storm::RationalNumber>::create;
 
-            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalNumber>> create(storm::dd::Bdd<DdType> 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;
+            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalNumber>> create() const override;
         };
         
         template<storm::dd::DdType DdType>
@@ -135,7 +147,7 @@ namespace storm {
         public:
             using SymbolicLinearEquationSolverFactory<DdType, storm::RationalFunction>::create;
 
-            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalFunction>> create(storm::dd::Bdd<DdType> 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;
+            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalFunction>> create() const override;
         };
         
     } // namespace solver
diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
index 490f341e3..4cb05ee90 100644
--- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
@@ -78,12 +78,12 @@ namespace storm {
         }
 
         template<storm::dd::DdType DdType, typename ValueType>
-        SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings) : settings(settings), requirementsChecked(false) {
+        SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings) : SymbolicEquationSolver<DdType, ValueType>(), settings(settings), requirementsChecked(false) {
             // Intentionally left empty.
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, storm::dd::Bdd<DdType> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::unique_ptr<SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory, SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings) : A(A), allRows(allRows), illegalMask(illegalMask), illegalMaskAdd(illegalMask.ite(A.getDdManager().getConstant(storm::utility::infinity<ValueType>()), A.getDdManager().template getAddZero<ValueType>())), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), choiceVariables(choiceVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), settings(settings), requirementsChecked(false) {
+        SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, storm::dd::Bdd<DdType> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::unique_ptr<SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory, SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings) : SymbolicEquationSolver<DdType, ValueType>(allRows), A(A), illegalMask(illegalMask), illegalMaskAdd(illegalMask.ite(A.getDdManager().getConstant(storm::utility::infinity<ValueType>()), A.getDdManager().template getAddZero<ValueType>())), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), choiceVariables(choiceVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), settings(settings), requirementsChecked(false) {
             // Intentionally left empty.
         }
         
@@ -104,7 +104,7 @@ namespace storm {
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        typename SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::ValueIterationResult SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion) const {
+        typename SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::ValueIterationResult SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const {
 
             // Set up local variables.
             storm::dd::Add<DdType, ValueType> localX = x;
@@ -112,7 +112,7 @@ namespace storm {
             
             // Value iteration loop.
             SolverStatus status = SolverStatus::InProgress;
-            while (status == SolverStatus::InProgress && iterations < this->settings.getMaximalNumberOfIterations()) {
+            while (status == SolverStatus::InProgress && iterations < maximalIterations) {
                 // Compute tmp = A * x + b
                 storm::dd::Add<DdType, ValueType> localXAsColumn = localX.swapVariables(this->rowColumnMetaVariablePairs);
                 storm::dd::Add<DdType, ValueType> tmp = this->A.multiplyMatrix(localXAsColumn, this->columnMetaVariables);
@@ -178,7 +178,7 @@ namespace storm {
         
         template<storm::dd::DdType DdType, typename ValueType>
         template<typename RationalType, typename ImpreciseType>
-        storm::dd::Add<DdType, RationalType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalX, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
+        storm::dd::Add<DdType, RationalType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
             
             // Storage for the rational sharpened vector.
             storm::dd::Add<DdType, RationalType> sharpenedX;
@@ -189,7 +189,7 @@ namespace storm {
             ValueType precision = this->getSettings().getPrecision();
             SolverStatus status = SolverStatus::InProgress;
             while (status == SolverStatus::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) {
-                typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion());
+                typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion(), this->settings.getMaximalNumberOfIterations());
                 
                 ++valueIterationInvocations;
                 STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << " completed in " << viResult.iterations << " iterations.");
@@ -226,19 +226,17 @@ namespace storm {
         template<storm::dd::DdType DdType, typename ValueType>
         template<typename ImpreciseType>
         typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
-            
-            return solveEquationsRationalSearchHelper<ValueType, ValueType>(dir, *this, *this, x, b, x, b);
+            return solveEquationsRationalSearchHelper<ValueType, ValueType>(dir, *this, *this, b, this->getLowerBounds(), b);
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
         template<typename ImpreciseType>
         typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
 
-            storm::dd::Add<DdType, storm::RationalNumber> rationalX = x.template toValueType<storm::RationalNumber>();
             storm::dd::Add<DdType, storm::RationalNumber> rationalB = b.template toValueType<storm::RationalNumber>();
             SymbolicMinMaxLinearEquationSolver<DdType, storm::RationalNumber> rationalSolver(this->A.template toValueType<storm::RationalNumber>(), this->allRows, this->illegalMask, this->rowMetaVariables, this->columnMetaVariables, this->choiceVariables, this->rowColumnMetaVariablePairs, std::make_unique<GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalNumber>>());
             
-            storm::dd::Add<DdType, storm::RationalNumber> rationalResult = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(dir, rationalSolver, *this, rationalX, rationalB, x, b);
+            storm::dd::Add<DdType, storm::RationalNumber> rationalResult = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(dir, rationalSolver, *this, rationalB, this->getLowerBounds(), b);
             return rationalResult.template toValueType<ValueType>();
         }
 
@@ -248,17 +246,18 @@ namespace storm {
             
             // First try to find a solution using the imprecise value type.
             storm::dd::Add<DdType, ValueType> rationalResult;
+            storm::dd::Add<DdType, ImpreciseType> impreciseX;
             try {
-                storm::dd::Add<DdType, ImpreciseType> impreciseX = x.template toValueType<ImpreciseType>();
+                impreciseX = this->getLowerBounds().template toValueType<ImpreciseType>();
                 storm::dd::Add<DdType, ImpreciseType> impreciseB = b.template toValueType<ImpreciseType>();
                 SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> impreciseSolver(this->A.template toValueType<ImpreciseType>(), this->allRows, this->illegalMask, this->rowMetaVariables, this->columnMetaVariables, this->choiceVariables, this->rowColumnMetaVariablePairs, std::make_unique<GeneralSymbolicLinearEquationSolverFactory<DdType, ImpreciseType>>());
                 
-                rationalResult = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(dir, *this, impreciseSolver, x, b, impreciseX, impreciseB);
+                rationalResult = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(dir, *this, impreciseSolver, b, impreciseX, impreciseB);
             } catch (storm::exceptions::PrecisionExceededException const& e) {
                 STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic.");
-                
+
                 // Fall back to precise value type if the precision of the imprecise value type was exceeded.
-                rationalResult = solveEquationsRationalSearchHelper<ValueType, ValueType>(dir, *this, *this, x, b, x, b);
+                rationalResult = solveEquationsRationalSearchHelper<ValueType, ValueType>(dir, *this, *this, b, impreciseX.template toValueType<ValueType>(), b);
             }
             return rationalResult.template toValueType<ValueType>();
         }
@@ -271,14 +270,16 @@ namespace storm {
         template<storm::dd::DdType DdType, typename ValueType>
         storm::dd::Add<DdType, ValueType>  SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
             // Set up the environment.
-            storm::dd::Add<DdType, ValueType> localX = x;
+            storm::dd::Add<DdType, ValueType> localX;
             
             // If we were given an initial scheduler, we take its solution as the starting point.
             if (this->hasInitialScheduler()) {
                 localX = solveEquationsWithScheduler(this->getInitialScheduler(), x, b);
+            } else {
+                localX = this->getLowerBounds();
             }
             
-            ValueIterationResult viResult = performValueIteration(dir, localX, b, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion());
+            ValueIterationResult viResult = performValueIteration(dir, localX, b, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion(), this->settings.getMaximalNumberOfIterations());
             
             if (viResult.status == SolverStatus::Converged) {
                 STORM_LOG_INFO("Iterative solver (value iteration) converged in " << viResult.iterations << " iterations.");
@@ -419,12 +420,14 @@ namespace storm {
                     }
                 }
             } else if (this->getSettings().getSolutionMethod() == SymbolicMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) {
+                requirements.requireLowerBounds();
                 if (equationSystemType == EquationSystemType::ReachabilityRewards) {
                     if (!direction || direction.get() == OptimizationDirection::Minimize) {
                         requirements.requireValidInitialScheduler();
                     }
                 }
             } else if (this->getSettings().getSolutionMethod() == SymbolicMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::RationalSearch) {
+                requirements.requireLowerBounds();
                 if (equationSystemType == EquationSystemType::ReachabilityRewards) {
                     if (!direction || direction.get() == OptimizationDirection::Minimize) {
                         requirements.requireNoEndComponents();
diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h
index 408faf21e..4e5fadbfd 100644
--- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h
+++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h
@@ -6,6 +6,7 @@
 #include <vector>
 #include <boost/optional.hpp>
 
+#include "storm/solver/SymbolicEquationSolver.h"
 #include "storm/solver/OptimizationDirection.h"
 #include "storm/solver/SymbolicLinearEquationSolver.h"
 #include "storm/solver/EquationSystemType.h"
@@ -58,7 +59,7 @@ namespace storm {
          * linear equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
          */
         template<storm::dd::DdType DdType, typename ValueType>
-        class SymbolicMinMaxLinearEquationSolver {
+        class SymbolicMinMaxLinearEquationSolver : public SymbolicEquationSolver<DdType, ValueType> {
         public:
             SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings = SymbolicMinMaxLinearEquationSolverSettings<ValueType>());
             
@@ -140,7 +141,7 @@ namespace storm {
             bool isRequirementsCheckedSet() const;
 
             /*!
-             * Determines whether the given vector x satisfies x = Ax + b.
+             * Determines whether the given vector x satisfies x = min/max Ax + b.
              */
             bool isSolution(OptimizationDirection dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
 
@@ -156,7 +157,7 @@ namespace storm {
             static storm::dd::Add<DdType, RationalType> sharpen(OptimizationDirection dir, uint64_t precision, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, RationalType> const& rationalB, bool& isSolution);
             
             template<typename RationalType, typename ImpreciseType>
-            storm::dd::Add<DdType, RationalType> solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalX, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const;
+            storm::dd::Add<DdType, RationalType> solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const;
             template<typename ImpreciseType>
             typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
             template<typename ImpreciseType>
@@ -177,15 +178,12 @@ namespace storm {
                 storm::dd::Add<DdType, ValueType> values;
             };
             
-            ValueIterationResult performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion) const;
+            ValueIterationResult performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const;
             
         protected:
             // The matrix defining the coefficients of the linear equation system.
             storm::dd::Add<DdType, ValueType> A;
             
-            // A BDD characterizing all rows of the equation system.
-            storm::dd::Bdd<DdType> allRows;
-            
             // A BDD characterizing the illegal choices.
             storm::dd::Bdd<DdType> illegalMask;
             
diff --git a/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp b/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
index 89f75bd41..5f3613100 100644
--- a/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
+++ b/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
@@ -9,6 +9,11 @@
 #include "storm/settings/SettingsManager.h"
 #include "storm/settings/modules/NativeEquationSolverSettings.h"
 
+#include "storm/utility/KwekMehlhorn.h"
+#include "storm/utility/macros.h"
+#include "storm/exceptions/NotSupportedException.h"
+#include "storm/exceptions/PrecisionExceededException.h"
+
 namespace storm {
     namespace solver {
         
@@ -18,6 +23,27 @@ namespace storm {
             storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>();
             
             // Get appropriate settings.
+            switch (settings.getLinearEquationSystemMethod()) {
+                case storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Power:
+                    this->method = SolutionMethod::Power;
+                    break;
+                case storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::RationalSearch:
+                    this->method = SolutionMethod::RationalSearch;
+                    break;
+                case storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Jacobi:
+                default:
+                    this->method = SolutionMethod::Jacobi;
+                    break;
+            }
+
+            // Adjust the method if none was specified and we are using rational numbers.
+            if (std::is_same<ValueType, storm::RationalNumber>::value) {
+                if (settings.isLinearEquationSystemTechniqueSetFromDefaultValue() && this->method != SolutionMethod::RationalSearch) {
+                    STORM_LOG_INFO("Selecting the rational search as the solution technique to guarantee exact results. If you want to override this, please explicitly specify a different method.");
+                    this->method = SolutionMethod::RationalSearch;
+                }
+            }
+            
             maximalNumberOfIterations = settings.getMaximalIterationCount();
             precision = storm::utility::convertNumber<ValueType>(settings.getPrecision());
             relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
@@ -38,6 +64,11 @@ namespace storm {
             this->relative = value;
         }
         
+        template<typename ValueType>
+        void SymbolicNativeLinearEquationSolverSettings<ValueType>::setSolutionMethod(SolutionMethod const& method) {
+            this->method = method;
+        }
+        
         template<typename ValueType>
         ValueType SymbolicNativeLinearEquationSolverSettings<ValueType>::getPrecision() const {
             return precision;
@@ -52,6 +83,16 @@ namespace storm {
         bool SymbolicNativeLinearEquationSolverSettings<ValueType>::getRelativeTerminationCriterion() const {
             return relative;
         }
+
+        template<typename ValueType>
+        typename SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod SymbolicNativeLinearEquationSolverSettings<ValueType>::getSolutionMethod() const {
+            return this->method;
+        }
+
+        template<storm::dd::DdType DdType, typename ValueType>
+        SymbolicNativeLinearEquationSolver<DdType, ValueType>::SymbolicNativeLinearEquationSolver(SymbolicNativeLinearEquationSolverSettings<ValueType> const& settings) : SymbolicLinearEquationSolver<DdType, ValueType>(), settings(settings) {
+            // Intentionally left empty.
+        }
         
         template<storm::dd::DdType DdType, typename ValueType>
         SymbolicNativeLinearEquationSolver<DdType, ValueType>::SymbolicNativeLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> 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, SymbolicNativeLinearEquationSolverSettings<ValueType> const& settings) : SymbolicNativeLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, settings) {
@@ -65,6 +106,18 @@ namespace storm {
         
         template<storm::dd::DdType DdType, typename ValueType>
         storm::dd::Add<DdType, ValueType>  SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
+                return solveEquationsJacobi(x, b);
+            } else if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Power) {
+                return solveEquationsPower(x, b);
+            } else if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::RationalSearch) {
+                return solveEquationsRationalSearch(x, b);
+            }
+            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The selected solution technique is not supported.");
+        }
+
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsJacobi(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
             storm::dd::DdManager<DdType>& manager = this->getDdManager();
             
             // Start by computing the Jacobi decomposition of the matrix A.
@@ -74,7 +127,7 @@ namespace storm {
             storm::dd::Add<DdType, ValueType> lu = diagonal.ite(manager.template getAddZero<ValueType>(), this->A);
             storm::dd::Add<DdType, ValueType> diagonalAdd = diagonal.template toAdd<ValueType>();
             storm::dd::Add<DdType, ValueType> diag = diagonalAdd.multiplyMatrix(this->A, this->columnMetaVariables);
-
+            
             storm::dd::Add<DdType, ValueType> scaledLu = lu / diag;
             storm::dd::Add<DdType, ValueType> scaledB = b / diag;
             
@@ -97,22 +150,197 @@ namespace storm {
             }
             
             if (converged) {
-                STORM_LOG_INFO("Iterative solver converged in " << iterationCount << " iterations.");
+                STORM_LOG_INFO("Iterative solver (jacobi) converged in " << iterationCount << " iterations.");
             } else {
-                STORM_LOG_WARN("Iterative solver did not converge in " << iterationCount << " iterations.");
+                STORM_LOG_WARN("Iterative solver (jacobi) did not converge in " << iterationCount << " iterations.");
             }
             
             return xCopy;
         }
         
+        template<storm::dd::DdType DdType, typename ValueType>
+        typename SymbolicNativeLinearEquationSolver<DdType, ValueType>::PowerIterationResult SymbolicNativeLinearEquationSolver<DdType, ValueType>::performPowerIteration(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const {
+            
+            // Set up additional environment variables.
+            storm::dd::Add<DdType, ValueType> currentX = x;
+            uint_fast64_t iterations = 0;
+            SolverStatus status = SolverStatus::InProgress;
+            
+            while (status == SolverStatus::InProgress && iterations < maximalIterations) {
+                storm::dd::Add<DdType, ValueType> currentXAsColumn = currentX.swapVariables(this->rowColumnMetaVariablePairs);
+                storm::dd::Add<DdType, ValueType> tmp = this->A.multiplyMatrix(currentXAsColumn, this->columnMetaVariables) + b;
+                
+                // Now check if the process already converged within our precision.
+                if (tmp.equalModuloPrecision(currentX, precision, relativeTerminationCriterion)) {
+                    status = SolverStatus::Converged;
+                }
+                
+                // Set up next iteration.
+                ++iterations;
+                currentX = tmp;
+            }
+
+            return PowerIterationResult(status, iterations, currentX);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsPower(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            PowerIterationResult result = performPowerIteration(this->getLowerBounds(), b, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion(), this->getSettings().getMaximalNumberOfIterations());
+            
+            if (result.status == SolverStatus::Converged) {
+                STORM_LOG_INFO("Iterative solver (power iteration) converged in " << result.iterations << " iterations.");
+            } else {
+                STORM_LOG_WARN("Iterative solver (power iteration) did not converge in " << result.iterations << " iterations.");
+            }
+            
+            return result.values;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        bool SymbolicNativeLinearEquationSolver<DdType, ValueType>::isSolutionFixedPoint(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            storm::dd::Add<DdType, ValueType> xAsColumn = x.swapVariables(this->rowColumnMetaVariablePairs);
+            storm::dd::Add<DdType, ValueType> tmp = this->A.multiplyMatrix(xAsColumn, this->columnMetaVariables);
+            tmp += b;
+
+            return x == tmp;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        template<typename RationalType, typename ImpreciseType>
+        storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::sharpen(uint64_t precision, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, RationalType> const& rationalB, bool& isSolution) {
+            
+            storm::dd::Add<DdType, RationalType> sharpenedX;
+            
+            for (uint64_t p = 1; p < precision; ++p) {
+                sharpenedX = x.sharpenKwekMehlhorn(precision);
+                isSolution = rationalSolver.isSolutionFixedPoint(sharpenedX, rationalB);
+                
+                if (isSolution) {
+                    break;
+                }
+            }
+            
+            return sharpenedX;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        template<typename RationalType, typename ImpreciseType>
+        storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
+            
+            // Storage for the rational sharpened vector.
+            storm::dd::Add<DdType, RationalType> sharpenedX;
+            
+            // The actual rational search.
+            uint64_t overallIterations = 0;
+            uint64_t powerIterationInvocations = 0;
+            ValueType precision = this->getSettings().getPrecision();
+            SolverStatus status = SolverStatus::InProgress;
+            while (status == SolverStatus::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) {
+                typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion(), this->settings.getMaximalNumberOfIterations());
+                
+                ++powerIterationInvocations;
+                STORM_LOG_TRACE("Completed " << powerIterationInvocations << " power iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations.");
+                
+                // Count the iterations.
+                overallIterations += result.iterations;
+                
+                // Compute maximal precision until which to sharpen.
+                uint64_t p = storm::utility::convertNumber<uint64_t>(storm::utility::ceil(storm::utility::log10<ValueType>(storm::utility::one<ValueType>() / precision)));
+                
+                bool isSolution = false;
+                sharpenedX = sharpen<RationalType, ImpreciseType>(p, rationalSolver, result.values, rationalB, isSolution);
+                
+                if (isSolution) {
+                    status = SolverStatus::Converged;
+                } else {
+                    precision = precision / 100;
+                }
+            }
+            
+            if (status == SolverStatus::InProgress) {
+                status = SolverStatus::MaximalIterationsExceeded;
+            }
+            
+            if (status == SolverStatus::Converged) {
+                STORM_LOG_INFO("Iterative solver (rational search) converged in " << overallIterations << " iterations.");
+            } else {
+                STORM_LOG_WARN("Iterative solver (rational search) did not converge in " << overallIterations << " iterations.");
+            }
+            
+            return sharpenedX;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        template<typename ImpreciseType>
+        typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            return solveEquationsRationalSearchHelper<ValueType, ValueType>(*this, *this, b, this->getLowerBounds(), b);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        template<typename ImpreciseType>
+        typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            
+            storm::dd::Add<DdType, storm::RationalNumber> rationalB = b.template toValueType<storm::RationalNumber>();
+            SymbolicNativeLinearEquationSolver<DdType, storm::RationalNumber> rationalSolver(this->A.template toValueType<storm::RationalNumber>(), this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs);
+            
+            storm::dd::Add<DdType, storm::RationalNumber> rationalResult = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(rationalSolver, *this, rationalB, this->getLowerBounds(), b);
+            return rationalResult.template toValueType<ValueType>();
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        template<typename ImpreciseType>
+        typename std::enable_if<!std::is_same<ValueType, ImpreciseType>::value, storm::dd::Add<DdType, ValueType>>::type SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            
+            // First try to find a solution using the imprecise value type.
+            storm::dd::Add<DdType, ValueType> rationalResult;
+            storm::dd::Add<DdType, ImpreciseType> impreciseX;
+            try {
+                impreciseX = this->getLowerBounds().template toValueType<ImpreciseType>();
+                storm::dd::Add<DdType, ImpreciseType> impreciseB = b.template toValueType<ImpreciseType>();
+                SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> impreciseSolver(this->A.template toValueType<ImpreciseType>(), this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs);
+                
+                rationalResult = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(*this, impreciseSolver, b, impreciseX, impreciseB);
+            } catch (storm::exceptions::PrecisionExceededException const& e) {
+                STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic.");
+                
+                // Fall back to precise value type if the precision of the imprecise value type was exceeded.
+                rationalResult = solveEquationsRationalSearchHelper<ValueType, ValueType>(*this, *this, b, impreciseX.template toValueType<ValueType>(), b);
+            }
+            return rationalResult.template toValueType<ValueType>();
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearch(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
+            return solveEquationsRationalSearchHelper<double>(x, b);
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        LinearEquationSolverProblemFormat SymbolicNativeLinearEquationSolver<DdType, ValueType>::getEquationProblemFormat() const {
+            if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
+                return LinearEquationSolverProblemFormat::EquationSystem;
+            }
+            return LinearEquationSolverProblemFormat::FixedPointSystem;
+        }
+        
+        template<storm::dd::DdType DdType, typename ValueType>
+        LinearEquationSolverRequirements SymbolicNativeLinearEquationSolver<DdType, ValueType>::getRequirements() const {
+            LinearEquationSolverRequirements requirements;
+            
+            if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Power || this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::RationalSearch) {
+                requirements.requireLowerBounds();
+            }
+            
+            return requirements;
+        }
+        
         template<storm::dd::DdType DdType, typename ValueType>
         SymbolicNativeLinearEquationSolverSettings<ValueType> const& SymbolicNativeLinearEquationSolver<DdType, ValueType>::getSettings() const {
             return settings;
         }
         
         template<storm::dd::DdType DdType, typename ValueType>
-        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicNativeLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Bdd<DdType> 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::make_unique<SymbolicNativeLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, settings);
+        std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicNativeLinearEquationSolverFactory<DdType, ValueType>::create() const {
+            return std::make_unique<SymbolicNativeLinearEquationSolver<DdType, ValueType>>(settings);
         }
     
         template<storm::dd::DdType DdType, typename ValueType>
diff --git a/src/storm/solver/SymbolicNativeLinearEquationSolver.h b/src/storm/solver/SymbolicNativeLinearEquationSolver.h
index 23b166a7c..b04b22fc5 100644
--- a/src/storm/solver/SymbolicNativeLinearEquationSolver.h
+++ b/src/storm/solver/SymbolicNativeLinearEquationSolver.h
@@ -1,6 +1,9 @@
 #pragma once
 
 #include "storm/solver/SymbolicLinearEquationSolver.h"
+#include "storm/solver/SolverStatus.h"
+
+#include "storm/utility/NumberTraits.h"
 
 namespace storm {
     namespace solver {
@@ -8,17 +11,26 @@ namespace storm {
         template<typename ValueType>
         class SymbolicNativeLinearEquationSolverSettings {
         public:
+            enum class SolutionMethod {
+                Jacobi, Power, RationalSearch
+            };
+            
             SymbolicNativeLinearEquationSolverSettings();
             
             void setPrecision(ValueType precision);
             void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations);
             void setRelativeTerminationCriterion(bool value);
+            void setSolutionMethod(SolutionMethod const& method);
             
             ValueType getPrecision() const;
             uint64_t getMaximalNumberOfIterations() const;
             bool getRelativeTerminationCriterion() const;
+            SolutionMethod getSolutionMethod() const;
             
         private:
+            // The selected solution method.
+            SolutionMethod method;
+            
             // The required precision for the iterative methods.
             ValueType precision;
             
@@ -37,6 +49,13 @@ namespace storm {
         template<storm::dd::DdType DdType, typename ValueType = double>
         class SymbolicNativeLinearEquationSolver : public SymbolicLinearEquationSolver<DdType, ValueType> {
         public:
+            /*!
+             * Constructs a symbolic linear equation solver.
+             *
+             * @param settings The settings to use.
+             */
+            SymbolicNativeLinearEquationSolver(SymbolicNativeLinearEquationSolverSettings<ValueType> const& settings = SymbolicNativeLinearEquationSolverSettings<ValueType>());
+            
             /*!
              * Constructs a symbolic linear equation solver with the given meta variable sets and pairs.
              *
@@ -73,11 +92,50 @@ namespace storm {
              * @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A.
              * @return The solution of the equation system.
              */
-            virtual storm::dd::Add<DdType, ValueType> solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
-
+            virtual storm::dd::Add<DdType, ValueType> solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const override;
+            
+            virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override;
+            virtual LinearEquationSolverRequirements getRequirements() const override;
+            
             SymbolicNativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
             
         private:
+            storm::dd::Add<DdType, ValueType> solveEquationsJacobi(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
+            storm::dd::Add<DdType, ValueType> solveEquationsPower(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
+            storm::dd::Add<DdType, ValueType> solveEquationsRationalSearch(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
+            
+            /*!
+             * Determines whether the given vector x satisfies x = Ax + b.
+             */
+            bool isSolutionFixedPoint(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
+            
+            template<typename RationalType, typename ImpreciseType>
+            static storm::dd::Add<DdType, RationalType> sharpen(uint64_t precision, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, RationalType> const& rationalB, bool& isSolution);
+            
+            template<typename RationalType, typename ImpreciseType>
+            storm::dd::Add<DdType, RationalType> solveEquationsRationalSearchHelper(SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const;
+            template<typename ImpreciseType>
+            typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
+            template<typename ImpreciseType>
+            typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
+            template<typename ImpreciseType>
+            typename std::enable_if<!std::is_same<ValueType, ImpreciseType>::value, storm::dd::Add<DdType, ValueType>>::type solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
+
+            template<storm::dd::DdType DdTypePrime, typename ValueTypePrime>
+            friend class SymbolicNativeLinearEquationSolver;
+
+            struct PowerIterationResult {
+                PowerIterationResult(SolverStatus status, uint64_t iterations, storm::dd::Add<DdType, ValueType> const& values) : status(status), iterations(iterations), values(values) {
+                    // Intentionally left empty.
+                }
+                
+                SolverStatus status;
+                uint64_t iterations;
+                storm::dd::Add<DdType, ValueType> values;
+            };
+            
+            PowerIterationResult performPowerIteration(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const;
+
             // The settings to use.
             SymbolicNativeLinearEquationSolverSettings<ValueType> settings;
         };
@@ -87,7 +145,7 @@ namespace storm {
         public:
             using SymbolicLinearEquationSolverFactory<DdType, ValueType>::create;
             
-            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> 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 override;
+            virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create() const override;
             
             SymbolicNativeLinearEquationSolverSettings<ValueType>& getSettings();
             SymbolicNativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
diff --git a/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp b/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp
index 6385c1d93..273ecb8fc 100644
--- a/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp
+++ b/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp
@@ -984,7 +984,11 @@ namespace storm {
                     end = std::chrono::high_resolution_clock::now();
                     
                     // Check quotient matrix for sanity.
-                    STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one<ValueType>()).isZero(), "Illegal entries in quotient matrix.");
+                    if (std::is_same<ValueType, storm::RationalNumber>::value) {
+                        STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one<ValueType>()).isZero(), "Illegal entries in quotient matrix.");
+                    } else {
+                        STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one<ValueType>() + storm::utility::convertNumber<ValueType>(1e-6)).isZero(), "Illegal entries in quotient matrix.");
+                    }
                     STORM_LOG_ASSERT(quotientTransitionMatrix.sumAbstract(blockPrimeVariableSet).equalModuloPrecision(quotientTransitionMatrix.notZero().existsAbstract(blockPrimeVariableSet).template toAdd<ValueType>(), ValueType(1e-6)), "Illegal non-probabilistic matrix.");
                     
                     STORM_LOG_TRACE("Quotient transition matrix extracted in " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << "ms.");
diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp
index 7fb973653..77665f6e5 100644
--- a/src/storm/utility/constants.cpp
+++ b/src/storm/utility/constants.cpp
@@ -708,6 +708,11 @@ namespace storm {
             return number;
         }
 
+        template<>
+        RationalFunction convertNumber(std::string const& number){
+            return RationalFunction(convertNumber<RationalFunctionCoefficient>(number));
+        }
+
         template<>
         RationalFunction convertNumber(storm::storage::sparse::state_type const& number) {
             return RationalFunction(convertNumber<RationalFunctionCoefficient>(number));