diff --git a/src/storm/settings/modules/MultiplierSettings.cpp b/src/storm/settings/modules/MultiplierSettings.cpp
index d386bbf3a..7e8058c45 100644
--- a/src/storm/settings/modules/MultiplierSettings.cpp
+++ b/src/storm/settings/modules/MultiplierSettings.cpp
@@ -15,7 +15,7 @@ namespace storm {
             const std::string MultiplierSettings::multiplierTypeOptionName = "type";
 
             MultiplierSettings::MultiplierSettings() : ModuleSettings(moduleName) {
-                std::vector<std::string> multiplierTypes = {"native", "gmmxx"};
+                std::vector<std::string> multiplierTypes = {"native", "inplace", "gmmxx"};
                 this->addOption(storm::settings::OptionBuilder(moduleName, multiplierTypeOptionName, true, "Sets which type of multiplier is preferred.")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a multiplier.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(multiplierTypes)).setDefaultValueString("gmmxx").build()).build());
                 
@@ -25,6 +25,8 @@ namespace storm {
                 std::string type = this->getOption(multiplierTypeOptionName).getArgumentByName("name").getValueAsString();
                 if (type == "native") {
                     return storm::solver::MultiplierType::Native;
+                } else if (type == "inplace") {
+                    return storm::solver::MultiplierType::InPlace;
                 } else if (type == "gmmxx") {
                     return storm::solver::MultiplierType::Gmmxx;
                 }
diff --git a/src/storm/solver/InPlaceMultiplier.cpp b/src/storm/solver/InPlaceMultiplier.cpp
new file mode 100644
index 000000000..abfd4463a
--- /dev/null
+++ b/src/storm/solver/InPlaceMultiplier.cpp
@@ -0,0 +1,138 @@
+#include "storm/solver/InPlaceMultiplier.h"
+
+#include "storm-config.h"
+
+#include "storm/environment/solver/MultiplierEnvironment.h"
+#include "storm/settings/SettingsManager.h"
+#include "storm/settings/modules/CoreSettings.h"
+
+#include "storm/storage/SparseMatrix.h"
+
+#include "storm/adapters/RationalNumberAdapter.h"
+#include "storm/adapters/RationalFunctionAdapter.h"
+
+#include "storm/utility/macros.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<typename ValueType>
+        InPlaceMultiplier<ValueType>::InPlaceMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix) : Multiplier<ValueType>(matrix) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        bool InPlaceMultiplier<ValueType>::parallelize(Environment const& env) const {
+#ifdef STORM_HAVE_INTELTBB
+            return storm::settings::getModule<storm::settings::modules::CoreSettings>().isUseIntelTbbSet();
+#else
+            return false;
+#endif
+        }
+        
+        template<typename ValueType>
+        void InPlaceMultiplier<ValueType>::multiply(Environment const& env, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
+            std::vector<ValueType>* target = &result;
+            if (&x == &result) {
+                if (this->cachedVector) {
+                    this->cachedVector->resize(x.size());
+                } else {
+                    this->cachedVector = std::make_unique<std::vector<ValueType>>(x.size());
+                }
+                target = this->cachedVector.get();
+            }
+            if (parallelize(env)) {
+                multAddParallel(x, b, *target);
+            } else {
+                multAdd(x, b, *target);
+            }
+            if (&x == &result) {
+                std::swap(result, *this->cachedVector);
+            }
+        }
+        
+        template<typename ValueType>
+        void InPlaceMultiplier<ValueType>::multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
+            this->matrix.multiplyWithVectorBackward(x, x, b);
+        }
+        
+        template<typename ValueType>
+        void InPlaceMultiplier<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 {
+            std::vector<ValueType>* target = &result;
+            if (&x == &result) {
+                if (this->cachedVector) {
+                    this->cachedVector->resize(x.size());
+                } else {
+                    this->cachedVector = std::make_unique<std::vector<ValueType>>(x.size());
+                }
+                target = this->cachedVector.get();
+            }
+            if (parallelize(env)) {
+                multAddReduceParallel(dir, rowGroupIndices, x, b, *target, choices);
+            } else {
+                multAddReduce(dir, rowGroupIndices, x, b, *target, choices);
+            }
+            if (&x == &result) {
+                std::swap(result, *this->cachedVector);
+            }
+        }
+        
+        template<typename ValueType>
+        void InPlaceMultiplier<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) const {
+            this->matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices);
+        }
+        
+        template<typename ValueType>
+        void InPlaceMultiplier<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 InPlaceMultiplier<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)) {
+                val1 += entry.getValue() * x1[entry.getColumn()];
+                val2 += entry.getValue() * x2[entry.getColumn()];
+            }
+        }
+
+
+        template<typename ValueType>
+        void InPlaceMultiplier<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 InPlaceMultiplier<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);
+        }
+        
+        template<typename ValueType>
+        void InPlaceMultiplier<ValueType>::multAddParallel(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
+#ifdef STORM_HAVE_INTELTBB
+            this->matrix.multiplyWithVectorParallel(x, result, b);
+#else
+            STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version.");
+            multAdd(x, b, result);
+#endif
+        }
+                
+        template<typename ValueType>
+        void InPlaceMultiplier<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 {
+#ifdef STORM_HAVE_INTELTBB
+            this->matrix.multiplyAndReduceParallel(dir, rowGroupIndices, x, b, result, choices);
+#else
+            STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version.");
+            multAddReduce(dir, rowGroupIndices, x, b, result, choices);
+#endif
+        }
+
+        template class InPlaceMultiplier<double>;
+#ifdef STORM_HAVE_CARL
+        template class InPlaceMultiplier<storm::RationalNumber>;
+        template class InPlaceMultiplier<storm::RationalFunction>;
+#endif
+        
+    }
+}
diff --git a/src/storm/solver/InPlaceMultiplier.h b/src/storm/solver/InPlaceMultiplier.h
new file mode 100644
index 000000000..d02c74bfb
--- /dev/null
+++ b/src/storm/solver/InPlaceMultiplier.h
@@ -0,0 +1,41 @@
+#pragma once
+
+#include "storm/solver/Multiplier.h"
+
+#include "storm/solver/OptimizationDirection.h"
+
+namespace storm {
+    namespace storage {
+        template<typename ValueType>
+        class SparseMatrix;
+    }
+    
+    namespace solver {
+        
+        template<typename ValueType>
+        class InPlaceMultiplier : public Multiplier<ValueType> {
+        public:
+            InPlaceMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix);
+            
+            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) 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) 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 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;
+            
+        };
+        
+    }
+}
diff --git a/src/storm/solver/Multiplier.cpp b/src/storm/solver/Multiplier.cpp
index 26bb4d85b..a57633db6 100644
--- a/src/storm/solver/Multiplier.cpp
+++ b/src/storm/solver/Multiplier.cpp
@@ -11,6 +11,7 @@
 #include "storm/solver/SolverSelectionOptions.h"
 #include "storm/solver/NativeMultiplier.h"
 #include "storm/solver/GmmxxMultiplier.h"
+#include "storm/solver/InPlaceMultiplier.h"
 #include "storm/environment/solver/MultiplierEnvironment.h"
 
 namespace storm {
@@ -63,6 +64,8 @@ namespace storm {
                     return std::make_unique<GmmxxMultiplier<ValueType>>(matrix);
                 case MultiplierType::Native:
                     return std::make_unique<NativeMultiplier<ValueType>>(matrix);
+                case MultiplierType::InPlace:
+                    return std::make_unique<InPlaceMultiplier<ValueType>>(matrix);
             }
         }
         
diff --git a/src/storm/solver/Multiplier.h b/src/storm/solver/Multiplier.h
index 87d51b2f5..77dc831c3 100644
--- a/src/storm/solver/Multiplier.h
+++ b/src/storm/solver/Multiplier.h
@@ -116,17 +116,17 @@ namespace storm {
              * Multiplies the row with the given index with x and adds the result to the provided value
              * @param rowIndex The index of the considered row
              * @param x The input vector with which the row is multiplied
-             * @param value The multiplication result is added to this value.
+             * @param value The multiplication result is added to this value. It shall not reffer to a value in x or in the Matrix.
              */
             virtual void multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType& value) const = 0;
             
             /*!
              * Multiplies the row with the given index with x1 and x2 and adds the given offset o1 and o2, respectively
              * @param rowIndex The index of the considered row
-             * @param x1 The first input vector with which the row is multiplied
-             * @param val1 The first multiplication result is added to this value.
-             * @param x2 The second input vector with which the row is multiplied
-             * @param val2 The second multiplication result is added to this value.
+             * @param x1 The first input vector with which the row is multiplied.
+             * @param val1 The first multiplication result is added to this value. It shall not reffer to a value in x or in the Matrix.
+             * @param x2 The second input vector with which the row is multiplied.
+             * @param val2 The second multiplication result is added to this value. It shall not reffer to a value in x or in the Matrix.
              */
             virtual void multiplyRow2(uint64_t const& rowIndex, std::vector<ValueType> const& x1, ValueType& val1, std::vector<ValueType> const& x2, ValueType& val2) const;
             
diff --git a/src/storm/solver/NativeMultiplier.cpp b/src/storm/solver/NativeMultiplier.cpp
index 13b27c617..f2e74bb71 100644
--- a/src/storm/solver/NativeMultiplier.cpp
+++ b/src/storm/solver/NativeMultiplier.cpp
@@ -1,3 +1,4 @@
+#include <storm/exceptions/InvalidTypeException.h>
 #include "storm/solver/NativeMultiplier.h"
 
 #include "storm-config.h"
@@ -12,17 +13,40 @@
 #include "storm/adapters/RationalFunctionAdapter.h"
 
 #include "storm/utility/macros.h"
+#include "storm/utility/vector.h"
 
 namespace storm {
     namespace solver {
         
-        template<typename ValueType>
-        NativeMultiplier<ValueType>::NativeMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix) : Multiplier<ValueType>(matrix) {
+        template<typename ValueType, typename IndexType>
+        NativeMultiplier<ValueType, IndexType>::NativeMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix) : Multiplier<ValueType>(matrix), numRows(0) {
             // Intentionally left empty.
         }
         
-        template<typename ValueType>
-        bool NativeMultiplier<ValueType>::parallelize(Environment const& env) const {
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::initialize() const {
+            if (numRows == 0) {
+                numRows = this->matrix.getRowCount();
+                entries.clear();
+                columns.clear();
+                rowIndications.clear();
+                entries.reserve(this->matrix.getNonzeroEntryCount());
+                columns.reserve(this->matrix.getColumnCount());
+                rowIndications.reserve(numRows + 1);
+                
+                rowIndications.push_back(0);
+                for (IndexType r = 0; r < numRows; ++r) {
+                    for (auto const& entry : this->matrix.getRow(r)) {
+                        entries.push_back(entry.getValue());
+                        columns.push_back(entry.getColumn());
+                    }
+                    rowIndications.push_back(entries.size());
+                }
+            }
+        }
+        
+        template<typename ValueType, typename IndexType>
+        bool NativeMultiplier<ValueType, IndexType>::parallelize(Environment const& env) const {
 #ifdef STORM_HAVE_INTELTBB
             return storm::settings::getModule<storm::settings::modules::CoreSettings>().isUseIntelTbbSet();
 #else
@@ -30,8 +54,9 @@ namespace storm {
 #endif
         }
         
-        template<typename ValueType>
-        void NativeMultiplier<ValueType>::multiply(Environment const& env, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::multiply(Environment const& env, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
+            initialize();
             std::vector<ValueType>* target = &result;
             if (&x == &result) {
                 if (this->cachedVector) {
@@ -51,13 +76,32 @@ namespace storm {
             }
         }
         
-        template<typename ValueType>
-        void NativeMultiplier<ValueType>::multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
-            this->matrix.multiplyWithVectorBackward(x, x, b);
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::multiplyGaussSeidel(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
+            initialize();
+            // multiply the rows in backwards order
+            IndexType r = numRows;
+            if (b) {
+                std::vector<ValueType> const& bRef = *b;
+                while (r > 0) {
+                    --r;
+                    ValueType xr = bRef[r];
+                    this->multiplyRow(r, x, xr);
+                    x[r] = std::move(xr);
+                }
+            } else {
+                while (r > 0) {
+                    --r;
+                    ValueType xr = storm::utility::zero<ValueType>();
+                    this->multiplyRow(r, x, xr);
+                    x[r] = std::move(xr);
+                }
+            }
         }
         
-        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 {
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::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 {
+            initialize();
             std::vector<ValueType>* target = &result;
             if (&x == &result) {
                 if (this->cachedVector) {
@@ -77,39 +121,132 @@ namespace storm {
             }
         }
         
-        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) const {
-            this->matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices);
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::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) const {
+            initialize();
+            assert(rowGroupIndices.size() - 1 == x.size());
+            // multiply the rowgroups in backwards order
+            IndexType g = x.size();
+            if (choices) {
+                while (g > 0) {
+                    --g;
+                    x[g] = multiplyAndReduceRowGroup(dir, rowGroupIndices[g], rowGroupIndices[g + 1], x, b, &((*choices)[g]));
+                }
+            } else {
+                while (g > 0) {
+                    --g;
+                    x[g] = multiplyAndReduceRowGroup(dir, rowGroupIndices[g], rowGroupIndices[g + 1], x, b);
+                }
+            }
         }
         
-        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, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType& value) const {
+            initialize();
+            assert(rowIndex < numRows);
+            IndexType const& rowStart = rowIndications[rowIndex];
+            IndexType const& rowEnd = rowIndications[rowIndex + 1];
+            for (IndexType e = rowStart; e < rowEnd; ++e) {
+                value += entries[e] * x[columns[e]];
             }
         }
         
-        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)) {
-                val1 += entry.getValue() * x1[entry.getColumn()];
-                val2 += entry.getValue() * x2[entry.getColumn()];
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::multiplyRow2(uint64_t const& rowIndex, std::vector<ValueType> const& x1, ValueType& val1, std::vector<ValueType> const& x2, ValueType& val2) const {
+            initialize();
+            assert(rowIndex < numRows);
+            IndexType const& rowStart = rowIndications[rowIndex];
+            IndexType const& rowEnd = rowIndications[rowIndex + 1];
+            for (IndexType e = rowStart; e < rowEnd; ++e) {
+                val1 += entries[e] * x1[columns[e]];
+                val2 += entries[e] * x2[columns[e]];
             }
         }
+        
+        template<typename ValueType, typename IndexType>
+        ValueType NativeMultiplier<ValueType, IndexType>::multiplyAndReduceRowGroup(OptimizationDirection const& dir, IndexType const& groupStart, IndexType const& groupEnd, std::vector<ValueType> const& x, std::vector<ValueType> const* b, uint_fast64_t* choice) const {
+            // Compute value for first row
+            ValueType result = b ? (*b)[groupStart] : storm::utility::zero<ValueType>();
+            multiplyRow(groupStart, x, result);
+            if (choice) {
+                *choice = 0;
+                // Compute the value for the remaining rows. Keep track of the optimal choice
+                for (IndexType r = groupStart + 1; r < groupEnd; ++r) {
+                    ValueType rowVal = b ? (*b)[r] : storm::utility::zero<ValueType>();
+                    multiplyRow(r, x, rowVal);
+                    if (dir == OptimizationDirection::Minimize) {
+                        if (rowVal < result) {
+                            result = std::move(rowVal);
+                            *choice = r - groupStart;
+                        }
+                    } else {
+                        if (rowVal > result) {
+                            result = std::move(rowVal);
+                            *choice = r - groupStart;
+                        }
+                    }
+                }
+            } else {
+                // Compute the value for the remaining rows
+                for (IndexType r = groupStart + 1; r < groupEnd; ++r) {
+                    ValueType rowVal = b ? (*b)[r] : storm::utility::zero<ValueType>();
+                    multiplyRow(r, x, rowVal);
+                    if (dir == OptimizationDirection::Minimize) {
+                        if (rowVal < result) {
+                            result = std::move(rowVal);
+                        }
+                    } else {
+                        if (rowVal > result) {
+                            result = std::move(rowVal);
+                        }
+                    }
+                }
+            }
+            return result;
+        }
+        
+        template<>
+        storm::RationalFunction NativeMultiplier<storm::RationalFunction, uint32_t>::multiplyAndReduceRowGroup(OptimizationDirection const&, uint32_t const&, uint32_t const&, std::vector<storm::RationalFunction> const&, std::vector<storm::RationalFunction> const*, uint_fast64_t*) const {
+            STORM_LOG_THROW(false, storm::exceptions::InvalidTypeException, "This operation is not supported with the given value type.");
+            return storm::utility::zero<storm::RationalFunction>();
+        }
 
-
-        template<typename ValueType>
-        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, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::multAdd(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
+            // multiply the rows sequentially (in forward order)
+            if (b) {
+                std::vector<ValueType> const& bRef = *b;
+                for (IndexType r = 0; r < numRows; ++r) {
+                    ValueType xr = bRef[r];
+                    this->multiplyRow(r, x, xr);
+                    result[r] = std::move(xr);
+                }
+            } else {
+                for (IndexType r = 0; r < numRows; ++r) {
+                    ValueType xr = storm::utility::zero<ValueType>();
+                    this->multiplyRow(r, x, xr);
+                    result[r] = std::move(xr);
+                }
+            }
         }
         
-        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);
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::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 {
+            assert(rowGroupIndices.size() - 1 == x.size());
+            // multiply the rowgroups in forward order
+            if (choices) {
+                for (IndexType g = 0, groupsEnd = x.size(); g < groupsEnd; ++g) {
+                    result[g] = multiplyAndReduceRowGroup(dir, rowGroupIndices[g], rowGroupIndices[g + 1], x, b, &((*choices)[g]));
+                }
+            } else {
+                for (IndexType g = 0, groupsEnd = x.size(); g < groupsEnd; ++g) {
+                    result[g] = multiplyAndReduceRowGroup(dir, rowGroupIndices[g], rowGroupIndices[g + 1], x, b);
+                }
+            }
         }
         
-        template<typename ValueType>
-        void NativeMultiplier<ValueType>::multAddParallel(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::multAddParallel(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
 #ifdef STORM_HAVE_INTELTBB
             this->matrix.multiplyWithVectorParallel(x, result, b);
 #else
@@ -118,8 +255,8 @@ namespace storm {
 #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 {
+        template<typename ValueType, typename IndexType>
+        void NativeMultiplier<ValueType, IndexType>::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 {
 #ifdef STORM_HAVE_INTELTBB
             this->matrix.multiplyAndReduceParallel(dir, rowGroupIndices, x, b, result, choices);
 #else
@@ -128,10 +265,10 @@ namespace storm {
 #endif
         }
 
-        template class NativeMultiplier<double>;
+        template class NativeMultiplier<double, uint32_t>;
 #ifdef STORM_HAVE_CARL
-        template class NativeMultiplier<storm::RationalNumber>;
-        template class NativeMultiplier<storm::RationalFunction>;
+        template class NativeMultiplier<storm::RationalNumber, uint32_t>;
+        template class NativeMultiplier<storm::RationalFunction, uint32_t>;
 #endif
         
     }
diff --git a/src/storm/solver/NativeMultiplier.h b/src/storm/solver/NativeMultiplier.h
index 73d1f84dc..d915b8fc7 100644
--- a/src/storm/solver/NativeMultiplier.h
+++ b/src/storm/solver/NativeMultiplier.h
@@ -12,7 +12,7 @@ namespace storm {
     
     namespace solver {
         
-        template<typename ValueType>
+        template<typename ValueType, typename IndexType = uint32_t>
         class NativeMultiplier : public Multiplier<ValueType> {
         public:
             NativeMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix);
@@ -28,6 +28,10 @@ namespace storm {
         private:
             bool parallelize(Environment const& env) const;
             
+            void initialize() const;
+            
+            ValueType multiplyAndReduceRowGroup(OptimizationDirection const& dir, IndexType const& groupStart, IndexType const& groupEnd, std::vector<ValueType> const& x, std::vector<ValueType> const* b, uint_fast64_t* choice = nullptr) 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;
@@ -35,6 +39,11 @@ namespace storm {
             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;
             
+            mutable std::vector<ValueType> entries;
+            mutable std::vector<IndexType> columns;
+            mutable std::vector<IndexType> rowIndications;
+            mutable IndexType numRows;
+            
         };
         
     }
diff --git a/src/storm/solver/SolverSelectionOptions.cpp b/src/storm/solver/SolverSelectionOptions.cpp
index d2148e9bc..bcbe9cc47 100644
--- a/src/storm/solver/SolverSelectionOptions.cpp
+++ b/src/storm/solver/SolverSelectionOptions.cpp
@@ -28,6 +28,8 @@ namespace storm {
             switch(t) {
                 case MultiplierType::Native:
                     return "Native";
+                case MultiplierType::InPlace:
+                    return "InPlace";
                 case MultiplierType::Gmmxx:
                     return "Gmmxx";
             }
diff --git a/src/storm/solver/SolverSelectionOptions.h b/src/storm/solver/SolverSelectionOptions.h
index 678bc1565..74d1dbc3c 100644
--- a/src/storm/solver/SolverSelectionOptions.h
+++ b/src/storm/solver/SolverSelectionOptions.h
@@ -7,7 +7,7 @@
 namespace storm {
     namespace solver {
         ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, LinearProgramming, Topological, RationalSearch, IntervalIteration, SoundValueIteration, TopologicalCuda)
-        ExtendEnumsWithSelectionField(MultiplierType, Native, Gmmxx)
+        ExtendEnumsWithSelectionField(MultiplierType, Native, InPlace, Gmmxx)
         ExtendEnumsWithSelectionField(GameMethod, PolicyIteration, ValueIteration)
         ExtendEnumsWithSelectionField(LraMethod, LinearProgramming, ValueIteration)