diff --git a/src/storm/solver/Multiplier.cpp b/src/storm/solver/Multiplier.cpp
index e48b29df2..b2b6dfd5e 100644
--- a/src/storm/solver/Multiplier.cpp
+++ b/src/storm/solver/Multiplier.cpp
@@ -10,6 +10,7 @@
 #include "storm/utility/macros.h"
 #include "storm/solver/SolverSelectionOptions.h"
 #include "storm/solver/NativeMultiplier.h"
+#include "storm/solver/GmmxxMultiplier.h"
 #include "storm/environment/solver/MultiplierEnvironment.h"
 
 namespace storm {
@@ -26,12 +27,12 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void Multiplier<ValueType>::multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices = nullptr) {
+        void Multiplier<ValueType>::multiplyAndReduce(Environment const& env, OptimizationDirection const& dir, std::vector<ValueType> const& x, std::vector<ValueType> const* b, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
             multiplyAndReduce(env, dir, this->matrix.getRowGroupIndices(), x, b, result, choices);
         }
 
         template<typename ValueType>
-        void Multiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices = nullptr) {
+        void Multiplier<ValueType>::multiplyAndReduceGaussSeidel(Environment const& env, OptimizationDirection const& dir, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint_fast64_t>* choices) const {
             multiplyAndReduceGaussSeidel(env, dir, this->matrix.getRowGroupIndices(), x, b, choices);
         }
     
@@ -43,9 +44,9 @@ namespace storm {
         }
     
         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 {
+        void Multiplier<ValueType>::repeatedMultiplyAndReduce(Environment const& env, OptimizationDirection const& dir, 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);
+                multiplyAndReduce(env, dir, x, b, x);
             }
         }
     
@@ -53,8 +54,7 @@ namespace storm {
         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");
+                    return std::make_unique<GmmxxMultiplier<ValueType>>(matrix);
                 case MultiplierType::Native:
                     return std::make_unique<NativeMultiplier<ValueType>>(matrix);
             }
diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp
index 47f07f6f8..c524513f2 100644
--- a/src/storm/solver/NativeLinearEquationSolver.cpp
+++ b/src/storm/solver/NativeLinearEquationSolver.cpp
@@ -301,11 +301,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        typename NativeLinearEquationSolver<ValueType>::PowerIterationResult NativeLinearEquationSolver<ValueType>::performPowerIteration(Environment const& env, std::vector<ValueType>*& currentX, std::vector<ValueType>*& newX, std::vector<ValueType> const& b, storm::solver::Multiplier<ValueType> const& multiplier, ValueType const& precision, bool relative, SolverGuarantee const& guarantee, uint64_t currentIterations, uint64_t maxIterations, storm::solver::MultiplicationStyle const& multiplicationStyle) const {
-
-            if (!this->multiplier) {
-                this->multiplier = storm::solver::MultiplierFactory<ValueType>().create(env, *A);
-            }
+        typename NativeLinearEquationSolver<ValueType>::PowerIterationResult NativeLinearEquationSolver<ValueType>::performPowerIteration(Environment const& env, std::vector<ValueType>*& currentX, std::vector<ValueType>*& newX, std::vector<ValueType> const& b, ValueType const& precision, bool relative, SolverGuarantee const& guarantee, uint64_t currentIterations, uint64_t maxIterations, storm::solver::MultiplicationStyle const& multiplicationStyle) const {
 
             bool useGaussSeidelMultiplication = multiplicationStyle == storm::solver::MultiplicationStyle::GaussSeidel;
             
@@ -317,9 +313,9 @@ namespace storm {
             while (!converged && !terminate && iterations < maxIterations) {
                 if (useGaussSeidelMultiplication) {
                     *newX = *currentX;
-                    multiplier.multiplyGaussSeidel(env, *newX, &b);
+                    this->multiplier->multiplyGaussSeidel(env, *newX, &b);
                 } else {
-                    multiplier.multiply(env, *currentX, &b, *newX);
+                    this->multiplier->multiply(env, *currentX, &b, *newX);
                 }
                 
                 // Now check for termination.
@@ -369,7 +365,7 @@ namespace storm {
             // Forward call to power iteration implementation.
             this->startMeasureProgress();
             ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision());
-            PowerIterationResult result = this->performPowerIteration(env, currentX, newX, b, *this->multiplier, precision, env.solver().native().getRelativeTerminationCriterion(), guarantee, 0, env.solver().native().getMaximalNumberOfIterations(), env.solver().native().getPowerMethodMultiplicationStyle());
+            PowerIterationResult result = this->performPowerIteration(env, currentX, newX, b, precision, env.solver().native().getRelativeTerminationCriterion(), guarantee, 0, env.solver().native().getMaximalNumberOfIterations(), env.solver().native().getPowerMethodMultiplicationStyle());
 
             // Swap the result in place.
             if (currentX == this->cachedRowVector.get()) {
@@ -596,7 +592,7 @@ namespace storm {
             
             void multiplyRow(uint64_t const& row, storm::storage::SparseMatrix<ValueType> const& A, storm::solver::Multiplier<ValueType> const& multiplier, ValueType const& bi, ValueType& xi, ValueType& yi) {
                 xi = multiplier.multiplyRow(row, x, bi);
-                yi = multiplier.multiplyRow(row, y, storm::utility::zero());
+                yi = multiplier.multiplyRow(row, y, storm::utility::zero<ValueType>());
                 /*
                 xi = bi;
                 yi = storm::utility::zero<ValueType>();
@@ -873,7 +869,7 @@ namespace storm {
             impreciseSolver.startMeasureProgress();
             while (status == SolverStatus::InProgress && overallIterations < maxIter) {
                 // Perform value iteration with the current precision.
-                typename NativeLinearEquationSolver<ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(currentX, newX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, SolverGuarantee::LessOrEqual, overallIterations, maxIter, multiplicationStyle);
+                typename NativeLinearEquationSolver<ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(env, currentX, newX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, SolverGuarantee::LessOrEqual, overallIterations, maxIter, multiplicationStyle);
                 
                 // At this point, the result of the imprecise value iteration is stored in the (imprecise) current x.
 
@@ -1143,7 +1139,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        LinearEquationSolverRequirements NativeLinearEquationSolver<ValueType>::getRequirements(Environment const& env, LinearEquationSolverTask const& task) const {
+        LinearEquationSolverRequirements NativeLinearEquationSolver<ValueType>::getRequirements(Environment const& env) const {
             LinearEquationSolverRequirements requirements;
             if (env.solver().native().isForceBoundsSet()) {
                 requirements.requireBounds();
diff --git a/src/storm/solver/NativeLinearEquationSolver.h b/src/storm/solver/NativeLinearEquationSolver.h
index 7640f19eb..46ba3628e 100644
--- a/src/storm/solver/NativeLinearEquationSolver.h
+++ b/src/storm/solver/NativeLinearEquationSolver.h
@@ -30,13 +30,6 @@ namespace storm {
             virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
             virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
             
-            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
-            virtual void multiplyAndReduce(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 bool supportsGaussSeidelMultiplication() const override;
-            virtual void multiplyGaussSeidel(std::vector<ValueType>& x, std::vector<ValueType> const* b) const override;
-            virtual void multiplyAndReduceGaussSeidel(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 ValueType multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x) const override;
-
             virtual LinearEquationSolverProblemFormat getEquationProblemFormat(storm::Environment const& env) const override;
             virtual LinearEquationSolverRequirements getRequirements(Environment const& env) const override;
 
diff --git a/src/storm/solver/NativeMultiplier.cpp b/src/storm/solver/NativeMultiplier.cpp
index 5fe7dde67..defe462f4 100644
--- a/src/storm/solver/NativeMultiplier.cpp
+++ b/src/storm/solver/NativeMultiplier.cpp
@@ -32,7 +32,6 @@ namespace storm {
         
         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 {
-            STORM_LOG_ASSERT(getMultiplicationStyle() == MultiplicationStyle::Regular, "Unexpected Multiplicationstyle.");
             std::vector<ValueType>* target = &result;
             if (&x == &result) {
                 if (this->cachedVector) {
@@ -53,13 +52,12 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void NativeMultiplier<ValueType>::multiplyGaussSeidel(Environment const& env, std::vector<ValueType> const& x, std::vector<ValueType> const* b) const {
+        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>
         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 {
-            STORM_LOG_ASSERT(getMultiplicationStyle() == MultiplicationStyle::Regular, "Unexpected Multiplicationstyle.");
             std::vector<ValueType>* target = &result;
             if (&x == &result) {
                 if (this->cachedVector) {
diff --git a/src/storm/solver/NativeMultiplier.h b/src/storm/solver/NativeMultiplier.h
index e91dc519e..e1cc4d01a 100644
--- a/src/storm/solver/NativeMultiplier.h
+++ b/src/storm/solver/NativeMultiplier.h
@@ -17,8 +17,6 @@ namespace storm {
         public:
             NativeMultiplier(storm::storage::SparseMatrix<ValueType> const& matrix);
             
-            virtual MultiplicationStyle getMultiplicationStyle() const override;
-
             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;
diff --git a/src/storm/solver/SolveGoal.h b/src/storm/solver/SolveGoal.h
index d93b86b91..9bd6f7e73 100644
--- a/src/storm/solver/SolveGoal.h
+++ b/src/storm/solver/SolveGoal.h
@@ -7,7 +7,6 @@
 #include "storm/solver/OptimizationDirection.h"
 #include "storm/logic/ComparisonType.h"
 #include "storm/storage/BitVector.h"
-#include "storm/solver/LinearEquationSolverTask.h"
 
 #include "storm/solver/LinearEquationSolver.h"
 #include "storm/solver/MinMaxLinearEquationSolver.h"
@@ -111,8 +110,8 @@ namespace storm {
         }
         
         template<typename ValueType, typename MatrixType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(Environment const& env, SolveGoal<ValueType>&& goal, storm::solver::LinearEquationSolverFactory<ValueType> const& factory, MatrixType&& matrix, storm::solver::LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) {
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = factory.create(env, std::forward<MatrixType>(matrix), task);
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> configureLinearEquationSolver(Environment const& env, SolveGoal<ValueType>&& goal, storm::solver::LinearEquationSolverFactory<ValueType> const& factory, MatrixType&& matrix) {
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = factory.create(env, std::forward<MatrixType>(matrix));
             if (goal.isBounded()) {
                 solver->setTerminationCondition(std::make_unique<TerminateIfFilteredExtremumExceedsThreshold<ValueType>>(goal.relevantValues(), goal.boundIsStrict(), goal.thresholdValue(), goal.minimize()));
             }
@@ -120,8 +119,8 @@ namespace storm {
         }
         
         template<typename MatrixType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> configureLinearEquationSolver(Environment const& env, SolveGoal<storm::RationalFunction>&& goal, storm::solver::LinearEquationSolverFactory<storm::RationalFunction> const& factory, MatrixType&& matrix, storm::solver::LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) {
-            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> solver = factory.create(env, std::forward<MatrixType>(matrix), task);
+        std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> configureLinearEquationSolver(Environment const& env, SolveGoal<storm::RationalFunction>&& goal, storm::solver::LinearEquationSolverFactory<storm::RationalFunction> const& factory, MatrixType&& matrix) {
+            std::unique_ptr<storm::solver::LinearEquationSolver<storm::RationalFunction>> solver = factory.create(env, std::forward<MatrixType>(matrix));
             return solver;
         }