diff --git a/src/storm/solver/Multiplier.cpp b/src/storm/solver/Multiplier.cpp
new file mode 100644
index 000000000..0375449bf
--- /dev/null
+++ b/src/storm/solver/Multiplier.cpp
@@ -0,0 +1,75 @@
+#include "storm/solver/Multiplier.h"
+
+#include "storm-config.h"
+
+#include "storm/storage/SparseMatrix.h"
+
+#include "storm/adapters/RationalNumberAdapter.h"
+#include "storm/adapters/RationalFunctionAdapter.h"
+
+#include "storm/utility/macros.h"
+#include "storm/solver/SolverSelectionOptions.h"
+#include "storm/solver/NativeMultiplier.h"
+#include "storm/environment/solver/MultiplierEnvironment.h"
+
+namespace storm {
+    namespace solver {
+        
+        template<typename ValueType>
+        Multiplier<ValueType>::Multiplier(storm::storage::SparseMatrix<ValueType> const& matrix) : matrix(matrix), allowGaussSeidelMultiplications(false) {
+            // Intentionally left empty.
+        }
+    
+        template<typename ValueType>
+        bool Multiplier<ValueType>::getAllowGaussSeidelMultiplications() const {
+            return allowGaussSeidelMultiplications;
+        }
+    
+        template<typename ValueType>
+        void Multiplier<ValueType>::setAllowGaussSeidelMultiplications(bool value) {
+            allowGaussSeidelMultiplications = value;
+        }
+        
+        template<typename ValueType>
+        void Multiplier<ValueType>::clearCache() const {
+            cachedVector.reset();
+        }
+    
+        template<typename ValueType>
+        void Multiplier<ValueType>::repeatedMultiply(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, uint64_t n) const {
+            for (uint64_t i = 0; i < n; ++i) {
+                multiply(env, x, b, x);
+            }
+        }
+    
+        template<typename ValueType>
+        void Multiplier<ValueType>::repeatedMultiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, uint64_t n) const {
+            for (uint64_t i = 0; i < n; ++i) {
+                multiplyAndReduce(env, dir, rowGroupIndices, x, b, x);
+            }
+        }
+    
+        template<typename ValueType>
+        std::unique_ptr<Multiplier<ValueType>> MultiplierFactory<ValueType>::create(Environment const& env, storm::storage::SparseMatrix<ValueType> const& matrix) {
+            switch (env.solver().multiplier().getType()) {
+                case MultiplierType::Gmmxx:
+                    //return std::make_unique<GmmxxMultiplier<ValueType>>(matrix);
+                STORM_PRINT_AND_LOG("gmm mult not yet supported");
+                case MultiplierType::Native:
+                    return std::make_unique<NativeMultiplier<ValueType>>(matrix);
+            }
+        }
+
+        
+        template class Multiplier<double>;
+        template class MultiplierFactory<double>;
+        
+#ifdef STORM_HAVE_CARL
+        template class Multiplier<storm::RationalNumber>;
+        template class MultiplierFactory<storm::RationalNumber>;
+        template class Multiplier<storm::RationalFunction>;
+        template class MultiplierFactory<storm::RationalFunction>;
+#endif
+        
+    }
+}
diff --git a/src/storm/solver/Multiplier.h b/src/storm/solver/Multiplier.h
new file mode 100644
index 000000000..432bce2ed
--- /dev/null
+++ b/src/storm/solver/Multiplier.h
@@ -0,0 +1,126 @@
+#pragma once
+
+#include "storm/solver/OptimizationDirection.h"
+#include "storm/solver/MultiplicationStyle.h"
+
+namespace storm {
+    
+    class Environment;
+    
+    namespace storage {
+        template<typename ValueType>
+        class SparseMatrix;
+    }
+    
+    namespace solver {
+        
+        template<typename ValueType>
+        class Multiplier {
+        public:
+            
+            Multiplier(storm::storage::SparseMatrix<ValueType> const& matrix);
+            
+            /*!
+             * Retrieves whether Gauss Seidel style multiplications are allowed.
+             */
+            bool getAllowGaussSeidelMultiplications() const;
+            
+            /*!
+             * Sets whether Gauss Seidel style multiplications are allowed.
+             */
+            void setAllowGaussSeidelMultiplications(bool value);
+            
+            /*!
+             * Returns the multiplication style performed by this multiplier
+             */
+            virtual MultiplicationStyle getMultiplicationStyle() const = 0;
+            
+            /*
+             * Clears the currently cached data of this multiplier in order to free some memory.
+             */
+            virtual void clearCache() const;
+            
+            /*!
+             * Performs a matrix-vector multiplication x' = A*x + b.
+             *
+             * @param x The input vector with which to multiply the matrix. Its length must be equal
+             * to the number of columns of A.
+             * @param b If non-null, this vector is added after the multiplication. If given, its length must be equal
+             * to the number of rows of A.
+             * @param result The target vector into which to write the multiplication result. Its length must be equal
+             * to the number of rows of A. Can be the same as the x vector.
+             */
+            virtual void multiply(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const = 0;
+            
+            /*!
+             * Performs a matrix-vector multiplication x' = A*x + b and then minimizes/maximizes over the row groups
+             * so that the resulting vector has the size of number of row groups of A.
+             *
+             * @param dir The direction for the reduction step.
+             * @param rowGroupIndices A vector storing the row groups over which to reduce.
+             * @param x The input vector with which to multiply the matrix. Its length must be equal
+             * to the number of columns of A.
+             * @param b If non-null, this vector is added after the multiplication. If given, its length must be equal
+             * to the number of rows of A.
+             * @param result The target vector into which to write the multiplication result. Its length must be equal
+             * to the number of rows of A. Can be the same as the x vector.
+             * @param choices If given, the choices made in the reduction process are written to this vector.
+             */
+            virtual void multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr) const = 0;
+            
+            /*!
+             * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After
+             * performing the necessary multiplications, the result is written to the input vector x. Note that the
+             * matrix A has to be given upon construction time of the solver object.
+             *
+             * @param x The initial vector with which to perform matrix-vector multiplication. Its length must be equal
+             * to the number of columns of A.
+             * @param b If non-null, this vector is added after each multiplication. If given, its length must be equal
+             * to the number of rows of A.
+             * @param n The number of times to perform the multiplication.
+             */
+            void repeatedMultiply(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, uint64_t n) const;
+            
+            /*!
+             * Performs repeated matrix-vector multiplication x' = A*x + b and then minimizes/maximizes over the row groups
+             * so that the resulting vector has the size of number of row groups of A.
+             *
+             * @param dir The direction for the reduction step.
+             * @param rowGroupIndices A vector storing the row groups over which to reduce.
+             * @param x The input vector with which to multiply the matrix. Its length must be equal
+             * to the number of columns of A.
+             * @param b If non-null, this vector is added after the multiplication. If given, its length must be equal
+             * to the number of rows of A.
+             * @param result The target vector into which to write the multiplication result. Its length must be equal
+             * to the number of rows of A.
+             * @param n The number of times to perform the multiplication.
+             */
+            void repeatedMultiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, uint64_t n) const;
+  
+            /*!
+             * Multiplies the row with the given index with x and adds the given offset
+             * @param rowIndex The index of the considered row
+             * @param x The input vector with which the row is multiplied
+             */
+            virtual ValueType multiplyRow(Environment const& env, uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType const& offset) const = 0;
+            
+        protected:
+            mutable std::unique_ptr<std::vector<ValueType>> cachedVector;
+            storm::storage::SparseMatrix<ValueType> const& matrix;
+        private:
+            bool allowGaussSeidelMultiplications;
+        };
+        
+        template<typename ValueType>
+        class MultiplierFactory {
+        public:
+            MultiplierFactory() = default;
+            ~MultiplierFactory() = default;
+
+            std::unique_ptr<Multiplier<ValueType>> create(Environment const& env, storm::storage::SparseMatrix<ValueType> const& matrix);
+            
+            
+        };
+        
+    }
+}
diff --git a/src/storm/solver/NativeMultiplier.cpp b/src/storm/solver/NativeMultiplier.cpp
index ddbeb235f..c4155fcbf 100644
--- a/src/storm/solver/NativeMultiplier.cpp
+++ b/src/storm/solver/NativeMultiplier.cpp
@@ -2,6 +2,10 @@
 
 #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"
@@ -13,89 +17,132 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        NativeMultiplier<ValueType>::NativeMultiplier() : storm::utility::VectorHelper<ValueType>() {
+        NativeMultiplier<ValueType>::NativeMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix) : Multiplier<ValueType>(matrix) {
             // Intentionally left empty.
         }
-
+        
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multAdd(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
-            std::vector<ValueType>* target = &result;
-            std::unique_ptr<std::vector<ValueType>> temporary;
-            if (&x == &result) {
-                STORM_LOG_WARN("Using temporary in 'multAdd'.");
-                temporary = std::make_unique<std::vector<ValueType>>(x.size());
-                target = temporary.get();
-            }
-            
-            if (this->parallelize()) {
-                multAddParallel(matrix, x, b, result);
+        bool NativeMultiplier<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>
+        MultiplicationStyle NativeMultiplier<ValueType>::getMultiplicationStyle() const {
+            if (this->getAllowGaussSeidelMultiplications()) {
+                return MultiplicationStyle::GaussSeidel;
             } else {
-                matrix.multiplyWithVector(x, result, b);
-            }
-            
-            if (target == temporary.get()) {
-                std::swap(result, *temporary);
+                return  MultiplicationStyle::Regular;
             }
         }
         
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multAddGaussSeidelBackward(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
-            matrix.multiplyWithVectorBackward(x, x, b);
+        void NativeMultiplier<ValueType>::multiply(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
+            if (getMultiplicationStyle() == MultiplicationStyle::GaussSeidel) {
+                multAddGaussSeidelBackward(x, b);
+                if (&x != &result) {
+                    result = x;
+                }
+            } else {
+                STORM_LOG_ASSERT(getMultiplicationStyle() == MultiplicationStyle::Regular, "Unexpected Multiplicationstyle.");
+                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 NativeMultiplier<ValueType>::multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices) const {
-            std::vector<ValueType>* target = &result;
-            std::unique_ptr<std::vector<ValueType>> temporary;
-            if (&x == &result) {
-                STORM_LOG_WARN("Using temporary in 'multAddReduce'.");
-                temporary = std::make_unique<std::vector<ValueType>>(x.size());
-                target = temporary.get();
-            }
-            
-            if (this->parallelize()) {
-                multAddReduceParallel(dir, rowGroupIndices, matrix, x, b, *target, choices);
+        void NativeMultiplier<ValueType>::multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+            if (getMultiplicationStyle() == MultiplicationStyle::GaussSeidel) {
+                multAddReduceGaussSeidelBackward(dir, rowGroupIndices, x, b, choices);
+                if (&x != &result) {
+                    result = x;
+                }
             } else {
-                matrix.multiplyAndReduce(dir, rowGroupIndices, x, b, *target, choices);
-            }
-            
-            if (target == temporary.get()) {
-                std::swap(result, *temporary);
+                STORM_LOG_ASSERT(getMultiplicationStyle() == MultiplicationStyle::Regular, "Unexpected Multiplicationstyle.");
+                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 NativeMultiplier<ValueType>::multAddReduceGaussSeidelBackward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint64_t>* choices) const {
-            matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices);
+        ValueType NativeMultiplier<ValueType>::multiplyRow(Environment const& env, uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType const& offset) const {
+            return this->matrix.multiplyRowWithVector(rowIndex, x);
+        }
+
+        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>
+        void NativeMultiplier<ValueType>::multAddGaussSeidelBackward(std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
+            this->matrix.multiplyWithVectorBackward(x, x, 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);
+        }
+        
+        template<typename ValueType>
+        void NativeMultiplier<ValueType>::multAddReduceGaussSeidelBackward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint64_t>* choices) const {
+            this->matrix.multiplyAndReduceBackward(dir, rowGroupIndices, x, b, x, choices);
         }
                 
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multAddParallel(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
+        void NativeMultiplier<ValueType>::multAddParallel(std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
 #ifdef STORM_HAVE_INTELTBB
             matrix.multiplyWithVectorParallel(x, result, b);
 #else
             STORM_LOG_WARN("Storm was built without support for Intel TBB, defaulting to sequential version.");
-            multAdd(matrix, x, b, result);
+            multAdd(x, b, result);
 #endif
         }
                 
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multAddReduceParallel(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, storm::storage::SparseMatrix<ValueType> const& matrix, 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) const {
 #ifdef STORM_HAVE_INTELTBB
             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, matrix, x, b, result, choices);
+            multAddReduce(dir, rowGroupIndices, x, b, result, choices);
 #endif
         }
 
-        template<typename ValueType>
-        ValueType NativeMultiplier<ValueType>::multiplyRow(storm::storage::SparseMatrix<ValueType> const& matrix, uint64_t const& rowIndex, std::vector<ValueType> const& x) const {
-            return matrix.multiplyRowWithVector(rowIndex, x);
-        }
-
         template class NativeMultiplier<double>;
-        
 #ifdef STORM_HAVE_CARL
         template class NativeMultiplier<storm::RationalNumber>;
         template class NativeMultiplier<storm::RationalFunction>;
diff --git a/src/storm/solver/NativeMultiplier.h b/src/storm/solver/NativeMultiplier.h
index ea8f543d4..b82598b0d 100644
--- a/src/storm/solver/NativeMultiplier.h
+++ b/src/storm/solver/NativeMultiplier.h
@@ -1,6 +1,6 @@
 #pragma once
 
-#include "storm/utility/VectorHelper.h"
+#include "storm/solver/Multiplier.h"
 
 #include "storm/solver/OptimizationDirection.h"
 
@@ -13,20 +13,28 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        class NativeMultiplier : public storm::utility::VectorHelper<ValueType> {
+        class NativeMultiplier : public Multiplier<ValueType> {
         public:
-            NativeMultiplier();
+            NativeMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix);
             
-            void multAdd(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const;
-            void multAddGaussSeidelBackward(storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType>& x, std::vector<ValueType> const* b) const;
+            virtual MultiplicationStyle getMultiplicationStyle() const override;
+
+            virtual void multiply(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
+            virtual void multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr) const override;
+            virtual ValueType multiplyRow(Environment const& env, uint64_t const& rowIndex, std::vector<ValueType> const& x, ValueType const& offset) 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 multAddGaussSeidelBackward(std::vector<ValueType>& x, std::vector<ValueType> const* b) const;
             
-            void multAddReduce(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint64_t>* choices = nullptr) const;
-            void multAddReduceGaussSeidelBackward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, storm::storage::SparseMatrix<ValueType> const& matrix, std::vector<ValueType>& x, std::vector<ValueType> const* b, 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) const;
+            void multAddReduceGaussSeidelBackward(storm::solver::OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint64_t>* choices = nullptr) const;
             
-            void multAddParallel(storm::storage::SparseMatrix<ValueType> const& matrix, 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, storm::storage::SparseMatrix<ValueType> const& matrix, 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;
             
-            ValueType multiplyRow(storm::storage::SparseMatrix<ValueType> const& matrix, uint64_t const& rowIndex, std::vector<ValueType> const& x) const;
         };
         
     }