diff --git a/src/storm/solver/NativeMultiplier.cpp b/src/storm/solver/NativeMultiplier.cpp
index 8b13d07cc..26dad0b22 100644
--- a/src/storm/solver/NativeMultiplier.cpp
+++ b/src/storm/solver/NativeMultiplier.cpp
@@ -60,9 +60,9 @@ namespace storm {
                 this->matrix.multiplyWithVectorForward(x, x, b);
             }
         }
-        
+
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+        void NativeMultiplier<ValueType>::multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
             std::vector<ValueType>* target = &result;
             if (&x == &result) {
                 if (this->cachedVector) {
@@ -73,31 +73,31 @@ namespace storm {
                 target = this->cachedVector.get();
             }
             if (parallelize(env)) {
-                multAddReduceParallel(dir, rowGroupIndices, x, b, *target, choices);
+                multAddReduceParallel(dir, rowGroupIndices, x, b, *target, choices, dirOverride);
             } else {
-                multAddReduce(dir, rowGroupIndices, x, b, *target, choices);
+                multAddReduce(dir, rowGroupIndices, x, b, *target, choices, dirOverride);
             }
             if (&x == &result) {
                 std::swap(result, *this->cachedVector);
             }
         }
-        
+
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices, bool backwards) const {
+        void NativeMultiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride, bool backwards) const {
             if (backwards) {
-                this->matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices);
+                this->matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices, dirOverride);
             } else {
-                this->matrix.multiplyAndReduceForward(dir, rowGroupIndices, x, b, x, choices);
+                this->matrix.multiplyAndReduceForward(dir, rowGroupIndices, x, b, x, choices, dirOverride);
             }
         }
-        
+
         template<typename ValueType>
         void NativeMultiplier<ValueType>::multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType& value) const {
             for (auto const& entry : this->matrix.getRow(rowIndex)) {
                 value += entry.getValue() * x[entry.getColumn()];
             }
         }
-        
+
         template<typename ValueType>
         void NativeMultiplier<ValueType>::multiplyRow2(uint64_t const& rowIndex, std::vector<ValueType> const& x1, ValueType& val1, std::vector<ValueType> const& x2, ValueType& val2) const {
             for (auto const& entry : this->matrix.getRow(rowIndex)) {
@@ -110,12 +110,12 @@ namespace storm {
         void NativeMultiplier<ValueType>::multAdd(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
             this->matrix.multiplyWithVector(x, result, b);
         }
-        
+
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices) const {
-            this->matrix.multiplyAndReduce(dir, rowGroupIndices, x, b, result, choices);
+        void NativeMultiplier<ValueType>::multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices, storm::storage::BitVector* dirOverride) const {
+            this->matrix.multiplyAndReduce(dir, rowGroupIndices, x, b, result, choices, dirOverride);
         }
-        
+
         template<typename ValueType>
         void NativeMultiplier<ValueType>::multAddParallel(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
 #ifdef STORM_HAVE_INTELTBB
@@ -125,11 +125,11 @@ namespace storm {
             multAdd(x, b, result);
 #endif
         }
-                
+
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices) const {
+        void NativeMultiplier<ValueType>::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices, storm::storage::BitVector* dirOverride) const {
 #ifdef STORM_HAVE_INTELTBB
-            this->matrix.multiplyAndReduceParallel(dir, rowGroupIndices, x, b, result, choices);
+            this->matrix.multiplyAndReduceParallel(dir, rowGroupIndices, x, b, result, choices, dirOverride);
 #else
             STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version.");
             multAddReduce(dir, rowGroupIndices, x, b, result, choices);
diff --git a/src/storm/solver/NativeMultiplier.h b/src/storm/solver/NativeMultiplier.h
index 92424f91c..62ae47d07 100644
--- a/src/storm/solver/NativeMultiplier.h
+++ b/src/storm/solver/NativeMultiplier.h
@@ -9,33 +9,33 @@ namespace storm {
         template<typename ValueType>
         class SparseMatrix;
     }
-    
+
     namespace solver {
-        
+
         template<typename ValueType>
         class NativeMultiplier : public Multiplier<ValueType> {
         public:
             NativeMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix);
             virtual ~NativeMultiplier() = default;
-            
+
             virtual void multiply(Environment const& env, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
             virtual void multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, bool backwards = true) const override;
-            virtual void multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr) const override;
-            virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr, bool backwards = true) const override;
+            virtual void multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr, storm::storage::BitVector* dirOverride = nullptr) const override;
+            virtual void multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr, storm::storage::BitVector* dirOverride = nullptr, bool backwards = true) const override;
             virtual void multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType& value) const override;
             virtual void multiplyRow2(uint64_t const& rowIndex, std::vector<ValueType> const& x1, ValueType& val1, std::vector<ValueType> const& x2, ValueType& val2) const override;
 
         private:
             bool parallelize(Environment const& env) const;
-            
+
             void multAdd(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const;
-            
-            void multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr) const;
-            
+
+            void multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr, storm::storage::BitVector* dirOverride = nullptr) const;
+
             void multAddParallel(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const;
-            void multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr) const;
-            
+            void multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr, storm::storage::BitVector* dirOverride = nullptr) const;
+
         };
-        
+
     }
 }
diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp
index d89f7a568..84b8fa306 100644
--- a/src/storm/storage/SparseMatrix.cpp
+++ b/src/storm/storage/SparseMatrix.cpp
@@ -1753,17 +1753,25 @@ namespace storm {
 #endif
 
         template<typename ValueType>
-        void SparseMatrix<ValueType>::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
-            if (dir == OptimizationDirection::Minimize) {
-                multiplyAndReduceForward<storm::utility::ElementLess<ValueType>>(rowGroupIndices, vector, summand, result, choices);
+        void SparseMatrix<ValueType>::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
+            if(dirOverride && !dirOverride->empty()) {
+                if (dir == OptimizationDirection::Minimize) {
+                    multiplyAndReduceForward<storm::utility::ElementLess<ValueType>, true>(rowGroupIndices, vector, summand, result, choices, dirOverride);
+                } else {
+                    multiplyAndReduceForward<storm::utility::ElementGreater<ValueType>, true>(rowGroupIndices, vector, summand, result, choices, dirOverride);
+                }
             } else {
-                multiplyAndReduceForward<storm::utility::ElementGreater<ValueType>>(rowGroupIndices, vector, summand, result, choices);
+                if (dir == OptimizationDirection::Minimize) {
+                    multiplyAndReduceForward<storm::utility::ElementLess<ValueType>, false>(rowGroupIndices, vector, summand, result, choices);
+                } else {
+                    multiplyAndReduceForward<storm::utility::ElementGreater<ValueType>, false>(rowGroupIndices, vector, summand, result, choices);
+                }
             }
         }
 
         template<typename ValueType>
-        template<typename Compare>
-        void SparseMatrix<ValueType>::multiplyAndReduceForward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+        template<typename Compare, bool dirOverridden>
+        void SparseMatrix<ValueType>::multiplyAndReduceForward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
             Compare compare;
             auto elementIt = this->begin();
             auto rowGroupIt = rowGroupIndices.begin();
@@ -1782,7 +1790,8 @@ namespace storm {
             uint64_t selectedChoice;
 
             uint64_t currentRow = 0;
-            for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++choiceIt, ++rowGroupIt) {
+            uint64_t currentRowGroup = 0;
+            for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++choiceIt, ++rowGroupIt, ++currentRowGroup) {
                 ValueType currentValue = storm::utility::zero<ValueType>();
 
                 // Only multiply and reduce if there is at least one row in the group.
@@ -1816,10 +1825,19 @@ namespace storm {
                             oldSelectedChoiceValue = newValue;
                         }
 
-                        if (compare(newValue, currentValue)) {
-                            currentValue = newValue;
-                            if (choices) {
-                                selectedChoice = currentRow - *rowGroupIt;
+                        if(dirOverridden) {
+                            if (dirOverride->get(currentRowGroup) ? compare(currentValue, newValue) : compare(newValue, currentValue)) {
+                                currentValue = newValue;
+                                if (choices) {
+                                    selectedChoice = currentRow - *rowGroupIt;
+                                }
+                            }
+                        } else {
+                            if (compare(newValue, currentValue)) {
+                                currentValue = newValue;
+                                if (choices) {
+                                    selectedChoice = currentRow - *rowGroupIt;
+                                }
                             }
                         }
                         if (summand) {
@@ -1829,8 +1847,14 @@ namespace storm {
 
                     // Finally write value to target vector.
                     *resultIt = currentValue;
-                    if (choices && compare(currentValue, oldSelectedChoiceValue)) {
-                        *choiceIt = selectedChoice;
+                    if(dirOverridden) {
+                        if (choices && dirOverride->get(currentRowGroup) ? compare(oldSelectedChoiceValue, currentValue) : compare(currentValue, oldSelectedChoiceValue)) {
+                            *choiceIt = selectedChoice;
+                        }
+                    } else {
+                        if (choices && compare(currentValue, oldSelectedChoiceValue)) {
+                            *choiceIt = selectedChoice;
+                        }
                     }
                 }
             }
@@ -1838,23 +1862,31 @@ namespace storm {
 
 #ifdef STORM_HAVE_CARL
         template<>
-        void SparseMatrix<storm::RationalFunction>::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction> const& vector, std::vector<storm::RationalFunction> const* b, std::vector<storm::RationalFunction>& result, std::vector<uint_fast64_t>* choices) const {
+        void SparseMatrix<storm::RationalFunction>::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction> const& vector, std::vector<storm::RationalFunction> const* b, std::vector<storm::RationalFunction>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
             STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported.");
         }
 #endif
 
         template<typename ValueType>
-        void SparseMatrix<ValueType>::multiplyAndReduceBackward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
-            if (dir == storm::OptimizationDirection::Minimize) {
-                multiplyAndReduceBackward<storm::utility::ElementLess<ValueType>>(rowGroupIndices, vector, summand, result, choices);
+        void SparseMatrix<ValueType>::multiplyAndReduceBackward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
+            if(dirOverride && !dirOverride->empty()) {
+                if (dir == storm::OptimizationDirection::Minimize) {
+                    multiplyAndReduceBackward<storm::utility::ElementLess<ValueType>, true>(rowGroupIndices, vector, summand, result, choices, dirOverride);
+                } else {
+                    multiplyAndReduceBackward<storm::utility::ElementGreater<ValueType>, true>(rowGroupIndices, vector, summand, result, choices, dirOverride);
+                }
             } else {
-                multiplyAndReduceBackward<storm::utility::ElementGreater<ValueType>>(rowGroupIndices, vector, summand, result, choices);
+                if (dir == storm::OptimizationDirection::Minimize) {
+                    multiplyAndReduceBackward<storm::utility::ElementLess<ValueType>, false>(rowGroupIndices, vector, summand, result, choices);
+                } else {
+                    multiplyAndReduceBackward<storm::utility::ElementGreater<ValueType>, false>(rowGroupIndices, vector, summand, result, choices);
+                }
             }
         }
 
         template<typename ValueType>
-        template<typename Compare>
-        void SparseMatrix<ValueType>::multiplyAndReduceBackward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+        template<typename Compare, bool directionOverridden>
+        void SparseMatrix<ValueType>::multiplyAndReduceBackward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
             Compare compare;
             auto elementIt = this->end() - 1;
             auto rowGroupIt = rowGroupIndices.end() - 2;
@@ -1873,7 +1905,8 @@ namespace storm {
             uint64_t selectedChoice;
 
             uint64_t currentRow = this->getRowCount() - 1;
-            for (auto resultIt = result.end() - 1, resultIte = result.begin() - 1; resultIt != resultIte; --resultIt, --choiceIt, --rowGroupIt) {
+            uint64_t currentRowGroup = this->getRowGroupCount() - 2;
+            for (auto resultIt = result.end() - 1, resultIte = result.begin() - 1; resultIt != resultIte; --resultIt, --choiceIt, --rowGroupIt, --currentRowGroup) {
                 ValueType currentValue = storm::utility::zero<ValueType>();
 
                 // Only multiply and reduce if there is at least one row in the group.
@@ -1905,18 +1938,33 @@ namespace storm {
                             oldSelectedChoiceValue = newValue;
                         }
 
-                        if (compare(newValue, currentValue)) {
-                            currentValue = newValue;
-                            if (choices) {
-                                selectedChoice = currentRow - *rowGroupIt;
+                        if(directionOverridden) {
+                            if (dirOverride->get(currentRowGroup) ? compare(currentValue, newValue) : compare(newValue, currentValue)) {
+                                currentValue = newValue;
+                                if (choices) {
+                                    selectedChoice = currentRow - *rowGroupIt;
+                                }
+                            }
+                        } else {
+                            if (compare(newValue, currentValue)) {
+                                currentValue = newValue;
+                                if (choices) {
+                                    selectedChoice = currentRow - *rowGroupIt;
+                                }
                             }
                         }
                     }
 
                     // Finally write value to target vector.
                     *resultIt = currentValue;
-                    if (choices && compare(currentValue, oldSelectedChoiceValue)) {
-                        *choiceIt = selectedChoice;
+                    if(directionOverridden) {
+                        if (choices && dirOverride->get(currentRowGroup) ? compare(oldSelectedChoiceValue, currentValue) : compare(currentValue, oldSelectedChoiceValue)) {
+                            *choiceIt = selectedChoice;
+                        }
+                    } else {
+                        if (choices && compare(currentValue, oldSelectedChoiceValue)) {
+                            *choiceIt = selectedChoice;
+                        }
                     }
                 }
             }
@@ -1924,20 +1972,20 @@ namespace storm {
 
 #ifdef STORM_HAVE_CARL
         template<>
-        void SparseMatrix<storm::RationalFunction>::multiplyAndReduceBackward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction> const& vector, std::vector<storm::RationalFunction> const* b, std::vector<storm::RationalFunction>& result, std::vector<uint_fast64_t>* choices) const {
+        void SparseMatrix<storm::RationalFunction>::multiplyAndReduceBackward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction> const& vector, std::vector<storm::RationalFunction> const* b, std::vector<storm::RationalFunction>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
             STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported.");
         }
 #endif
 
 #ifdef STORM_HAVE_INTELTBB
-        template <typename ValueType, typename Compare>
+        template <typename ValueType, typename Compare, bool directionOverridden>
         class TbbMultAddReduceFunctor {
         public:
             typedef typename storm::storage::SparseMatrix<ValueType>::index_type index_type;
             typedef typename storm::storage::SparseMatrix<ValueType>::value_type value_type;
             typedef typename storm::storage::SparseMatrix<ValueType>::const_iterator const_iterator;
 
-            TbbMultAddReduceFunctor(std::vector<uint64_t> const& rowGroupIndices, std::vector<MatrixEntry<index_type, value_type>> const& columnsAndEntries, std::vector<uint64_t> const& rowIndications, std::vector<ValueType> const& x, std::vector<ValueType>& result, std::vector<value_type> const* summand, std::vector<uint_fast64_t>* choices) : rowGroupIndices(rowGroupIndices), columnsAndEntries(columnsAndEntries), rowIndications(rowIndications), x(x), result(result), summand(summand), choices(choices) {
+            TbbMultAddReduceFunctor(std::vector<uint64_t> const& rowGroupIndices, std::vector<MatrixEntry<index_type, value_type>> const& columnsAndEntries, std::vector<uint64_t> const& rowIndications, std::vector<ValueType> const& x, std::vector<ValueType>& result, std::vector<value_type> const* summand, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride = nullptr) : rowGroupIndices(rowGroupIndices), columnsAndEntries(columnsAndEntries), rowIndications(rowIndications), x(x), result(result), summand(summand), choices(choices), dirOverride(dirOverride) {
                 // Intentionally left empty.
             }
 
@@ -1964,7 +2012,8 @@ namespace storm {
                 uint64_t selectedChoice;
 
                 uint64_t currentRow = *groupIt;
-                for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) {
+                uint64_t currentRowGroup = *groupIt;
+                for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt, ++currentRowGroup) {
                     ValueType currentValue = storm::utility::zero<ValueType>();
 
                     // Only multiply and reduce if there is at least one row in the group.
@@ -1998,18 +2047,33 @@ namespace storm {
                                 oldSelectedChoiceValue = newValue;
                             }
 
-                            if (compare(newValue, currentValue)) {
-                                currentValue = newValue;
-                                if (choices) {
-                                    selectedChoice = currentRow - *groupIt;
+                            if(directionOverridden) {
+                                if (dirOverride.get()->get(currentRowGroup) ? compare(currentValue, newValue) : compare(newValue, currentValue)) {
+                                    currentValue = newValue;
+                                    if (choices) {
+                                        selectedChoice = currentRow - *groupIt;
+                                    }
+                                }
+                            } else {
+                                if (compare(newValue, currentValue)) {
+                                    currentValue = newValue;
+                                    if (choices) {
+                                        selectedChoice = currentRow - *groupIt;
+                                    }
                                 }
                             }
                         }
 
                         // Finally write value to target vector.
                         *resultIt = currentValue;
-                        if (choices && compare(currentValue, oldSelectedChoiceValue)) {
-                            *choiceIt = selectedChoice;
+                        if(directionOverridden) {
+                            if (choices && dirOverride.get()->get(currentRowGroup) ? compare(oldSelectedChoiceValue, currentValue) : compare(currentValue, oldSelectedChoiceValue)) {
+                                *choiceIt = selectedChoice;
+                            }
+                        } else {
+                            if (choices && compare(currentValue, oldSelectedChoiceValue)) {
+                                *choiceIt = selectedChoice;
+                            }
                         }
                     }
                 }
@@ -2024,27 +2088,36 @@ namespace storm {
             std::vector<ValueType>& result;
             std::vector<value_type> const* summand;
             std::vector<uint_fast64_t>* choices;
+            boost::optional<storm::storage::BitVector*> dirOverride = boost::none;
         };
 
         template<typename ValueType>
-        void SparseMatrix<ValueType>::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
-            if (dir == storm::OptimizationDirection::Minimize) {
-                tbb::parallel_for(tbb::blocked_range<index_type>(0, rowGroupIndices.size() - 1, 100), TbbMultAddReduceFunctor<ValueType, storm::utility::ElementLess<ValueType>>(rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices));
+        void SparseMatrix<ValueType>::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
+            if(dirOverride && !dirOverride->empty()) {
+                if (dir == storm::OptimizationDirection::Minimize) {
+                    tbb::parallel_for(tbb::blocked_range<index_type>(0, rowGroupIndices.size() - 1, 100), TbbMultAddReduceFunctor<ValueType, storm::utility::ElementLess<ValueType>, true>(rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices, dirOverride));
+                } else {
+                    tbb::parallel_for(tbb::blocked_range<index_type>(0, rowGroupIndices.size() - 1, 100), TbbMultAddReduceFunctor<ValueType, storm::utility::ElementGreater<ValueType>, true>(rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices, dirOverride));
+                }
             } else {
-                tbb::parallel_for(tbb::blocked_range<index_type>(0, rowGroupIndices.size() - 1, 100), TbbMultAddReduceFunctor<ValueType, storm::utility::ElementGreater<ValueType>>(rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices));
+                if (dir == storm::OptimizationDirection::Minimize) {
+                    tbb::parallel_for(tbb::blocked_range<index_type>(0, rowGroupIndices.size() - 1, 100), TbbMultAddReduceFunctor<ValueType, storm::utility::ElementLess<ValueType>, false>(rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices));
+                } else {
+                    tbb::parallel_for(tbb::blocked_range<index_type>(0, rowGroupIndices.size() - 1, 100), TbbMultAddReduceFunctor<ValueType, storm::utility::ElementGreater<ValueType>, false>(rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices));
+                }
             }
         }
 
 #ifdef STORM_HAVE_CARL
         template<>
-        void SparseMatrix<storm::RationalFunction>::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction> const& vector, std::vector<storm::RationalFunction> const* summand, std::vector<storm::RationalFunction>& result, std::vector<uint_fast64_t>* choices) const {
+        void SparseMatrix<storm::RationalFunction>::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction> const& vector, std::vector<storm::RationalFunction> const* summand, std::vector<storm::RationalFunction>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
             STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported.");
         }
 #endif
 #endif
 
         template<typename ValueType>
-        void SparseMatrix<ValueType>::multiplyAndReduce(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+        void SparseMatrix<ValueType>::multiplyAndReduce(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride) const {
 
             // If the vector and the result are aliases, we need and temporary vector.
             std::vector<ValueType>* target;
diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h
index ba25a2819..8346e5554 100644
--- a/src/storm/storage/SparseMatrix.h
+++ b/src/storm/storage/SparseMatrix.h
@@ -907,19 +907,19 @@ namespace storm {
              * the 'new' choice has a value strictly better (wrt. to the optimization direction) value.
              * @return The resulting vector the content of the given result vector.
              */
-            void multiplyAndReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const;
-            
-            void multiplyAndReduceForward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const;
-            template<typename Compare>
-            void multiplyAndReduceForward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const;
-            
-            void multiplyAndReduceBackward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const;
-            template<typename Compare>
-            void multiplyAndReduceBackward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const;
+            void multiplyAndReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride = nullptr) const;
+
+            void multiplyAndReduceForward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride = nullptr) const;
+            template<typename Compare, bool directionOverridden>
+            void multiplyAndReduceForward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride = nullptr) const;
+
+            void multiplyAndReduceBackward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride = nullptr) const;
+            template<typename Compare, bool directionOverridden>
+            void multiplyAndReduceBackward(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride = nullptr) const;
 #ifdef STORM_HAVE_INTELTBB
-            void multiplyAndReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const;
-            template<typename Compare>
-            void multiplyAndReduceParallel(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const;
+            void multiplyAndReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride = nullptr) const;
+            template<typename Compare, bool directionOverridden>
+            void multiplyAndReduceParallel(std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices, storm::storage::BitVector* dirOverride = nullptr) const;
 #endif
 
             /*!