diff --git a/src/storm/solver/EigenLinearEquationSolver.cpp b/src/storm/solver/EigenLinearEquationSolver.cpp
index bcba06888..afbc88a3b 100644
--- a/src/storm/solver/EigenLinearEquationSolver.cpp
+++ b/src/storm/solver/EigenLinearEquationSolver.cpp
@@ -240,41 +240,6 @@ namespace storm {
             return true;
         }
         
-        template<typename ValueType>
-        void EigenLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
-            // Typedef the map-type so we don't have to spell it out.
-            typedef decltype(StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
-
-            auto eigenX = StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(x.data(), x.size());
-            auto eigenResult = StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(result.data(), result.size());
-
-            std::unique_ptr<MapType> eigenB;
-            if (b != nullptr) {
-                eigenB = std::make_unique<MapType>(StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(b->data(), b->size()));
-            }
-            
-            if (&x != &result) {
-                if (b != nullptr) {
-                    eigenResult.noalias() = *eigenA * eigenX + *eigenB;
-                } else {
-                    eigenResult.noalias() = *eigenA * eigenX;
-                }
-            } else {
-                if (b != nullptr) {
-                    eigenResult = *eigenA * eigenX + *eigenB;
-                } else {
-                    eigenResult = *eigenA * eigenX;
-                }
-            }
-        }
-        
-        template<typename ValueType>
-        ValueType EigenLinearEquationSolver<ValueType>::multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x) const {
-            auto eigenX = StormEigen::Matrix<ValueType, StormEigen::Dynamic, 1>::Map(x.data(), x.size());
-            return (eigenA->row(rowIndex) * eigenX)(0);
-        }
-
-        
         template<typename ValueType>
         LinearEquationSolverProblemFormat EigenLinearEquationSolver<ValueType>::getEquationProblemFormat(Environment const& env) const {
             return LinearEquationSolverProblemFormat::EquationSystem;
@@ -291,7 +256,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EigenLinearEquationSolverFactory<ValueType>::create(Environment const& env, LinearEquationSolverTask const& task) const {
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EigenLinearEquationSolverFactory<ValueType>::create(Environment const& env) const {
             return std::make_unique<storm::solver::EigenLinearEquationSolver<ValueType>>();
         }
         
diff --git a/src/storm/solver/EigenLinearEquationSolver.h b/src/storm/solver/EigenLinearEquationSolver.h
index 0bdeabba1..63768563a 100644
--- a/src/storm/solver/EigenLinearEquationSolver.h
+++ b/src/storm/solver/EigenLinearEquationSolver.h
@@ -20,10 +20,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 ValueType multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x) const override;
-
-
             virtual LinearEquationSolverProblemFormat getEquationProblemFormat(Environment const& env) const override;
 
         protected:
@@ -45,7 +41,7 @@ namespace storm {
         public:
             using LinearEquationSolverFactory<ValueType>::create;
 
-            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(Environment const& env, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(Environment const& env) const override;
 
             virtual std::unique_ptr<LinearEquationSolverFactory<ValueType>> clone() const override;
         };
diff --git a/src/storm/solver/EliminationLinearEquationSolver.cpp b/src/storm/solver/EliminationLinearEquationSolver.cpp
index 2e67119d6..2873c04c4 100644
--- a/src/storm/solver/EliminationLinearEquationSolver.cpp
+++ b/src/storm/solver/EliminationLinearEquationSolver.cpp
@@ -93,28 +93,6 @@ namespace storm {
             return true;
         }
         
-        template<typename ValueType>
-        void EliminationLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
-            if (&x != &result) {
-                A->multiplyWithVector(x, result);
-                if (b != nullptr) {
-                    storm::utility::vector::addVectors(result, *b, result);
-                }
-            } else {
-                // If the two vectors are aliases, we need to create a temporary.
-                std::vector<ValueType> tmp(result.size());
-                A->multiplyWithVector(x, tmp);
-                if (b != nullptr) {
-                    storm::utility::vector::addVectors(tmp, *b, result);
-                }
-            }
-        }
-
-        template<typename ValueType>
-        ValueType EliminationLinearEquationSolver<ValueType>::multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x) const {
-            return A->multiplyRowWithVector(rowIndex, x);
-        }
-        
         template<typename ValueType>
         LinearEquationSolverProblemFormat EliminationLinearEquationSolver<ValueType>::getEquationProblemFormat(Environment const& env) const {
             return LinearEquationSolverProblemFormat::FixedPointSystem;
@@ -131,7 +109,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EliminationLinearEquationSolverFactory<ValueType>::create(Environment const& env, LinearEquationSolverTask const& task) const {
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EliminationLinearEquationSolverFactory<ValueType>::create(Environment const& env) const {
             return std::make_unique<storm::solver::EliminationLinearEquationSolver<ValueType>>();
         }
         
diff --git a/src/storm/solver/EliminationLinearEquationSolver.h b/src/storm/solver/EliminationLinearEquationSolver.h
index 57da6bbaa..d3c3f5992 100644
--- a/src/storm/solver/EliminationLinearEquationSolver.h
+++ b/src/storm/solver/EliminationLinearEquationSolver.h
@@ -21,10 +21,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 ValueType multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x) const override;
-
-
             virtual LinearEquationSolverProblemFormat getEquationProblemFormat(Environment const& env) const override;
 
         protected:
@@ -48,7 +44,7 @@ namespace storm {
         public:
             using LinearEquationSolverFactory<ValueType>::create;
             
-            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(Environment const& env, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(Environment const& env) const override;
             
             virtual std::unique_ptr<LinearEquationSolverFactory<ValueType>> clone() const override;
 
diff --git a/src/storm/solver/GmmxxLinearEquationSolver.cpp b/src/storm/solver/GmmxxLinearEquationSolver.cpp
index f87a02a96..83d5b6ef6 100644
--- a/src/storm/solver/GmmxxLinearEquationSolver.cpp
+++ b/src/storm/solver/GmmxxLinearEquationSolver.cpp
@@ -5,7 +5,6 @@
 
 #include "storm/adapters/GmmxxAdapter.h"
 
-#include "storm/solver/GmmxxMultiplier.h"
 #include "storm/environment/solver/GmmxxSolverEnvironment.h"
 
 #include "storm/utility/vector.h"
@@ -116,41 +115,6 @@ namespace storm {
             return false;
         }
         
-        template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
-            multiplier.multAdd(*gmmxxA, x, b, result);
-            
-            if (!this->isCachingEnabled()) {
-                clearCache();
-            }
-        }
-        
-        template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::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) const {
-            multiplier.multAddReduce(dir, rowGroupIndices, *gmmxxA, x, b, result, choices);
-        }
-        
-        template<typename ValueType>
-        bool GmmxxLinearEquationSolver<ValueType>::supportsGaussSeidelMultiplication() const {
-            return true;
-        }
-        
-        template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::multiplyGaussSeidel(std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
-            multiplier.multAddGaussSeidelBackward(*gmmxxA, x, b);
-        }
-
-        template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::multiplyAndReduceGaussSeidel(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<uint64_t>* choices) const {
-            multiplier.multAddReduceGaussSeidel(dir, rowGroupIndices, *gmmxxA, x, b, choices);
-        }
-        
-        template<typename ValueType>
-        ValueType GmmxxLinearEquationSolver<ValueType>::multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x) const {
-            return multiplier.multiplyRow(*gmmxxA, rowIndex, x);
-        }
-
-        
         template<typename ValueType>
         LinearEquationSolverProblemFormat GmmxxLinearEquationSolver<ValueType>::getEquationProblemFormat(Environment const& env) const {
             return LinearEquationSolverProblemFormat::EquationSystem;
@@ -174,7 +138,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(Environment const& env, LinearEquationSolverTask const& task) const {
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(Environment const& env) const {
             return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>();
         }
         
diff --git a/src/storm/solver/GmmxxLinearEquationSolver.h b/src/storm/solver/GmmxxLinearEquationSolver.h
index 434d5688a..f288e9f58 100644
--- a/src/storm/solver/GmmxxLinearEquationSolver.h
+++ b/src/storm/solver/GmmxxLinearEquationSolver.h
@@ -5,10 +5,8 @@
 
 #include "storm/utility/gmm.h"
 
-#include "storm/solver/GmmxxMultiplier.h"
-
 #include "storm/solver/LinearEquationSolver.h"
-#include "SolverSelectionOptions.h"
+#include "storm/solver/SolverSelectionOptions.h"
 
 namespace storm {
     namespace solver {
@@ -27,13 +25,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(Environment const& env) const override;
             
             virtual void clearCache() const override;
@@ -51,9 +42,6 @@ namespace storm {
             // The matrix in gmm++ format.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxA;
             
-            // A multiplier object used to dispatch the multiplication calls.
-            GmmxxMultiplier<ValueType> multiplier;
-            
             // cached data obtained during solving
             mutable std::unique_ptr<gmm::ilu_precond<gmm::csr_matrix<ValueType>>> iluPreconditioner;
             mutable std::unique_ptr<gmm::diagonal_precond<gmm::csr_matrix<ValueType>>> diagonalPreconditioner;
@@ -64,7 +52,7 @@ namespace storm {
         public:
             using LinearEquationSolverFactory<ValueType>::create;
             
-            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(Environment const& env, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(Environment const& env) const override;
             
             virtual std::unique_ptr<LinearEquationSolverFactory<ValueType>> clone() const override;
 
diff --git a/src/storm/solver/LinearEquationSolver.cpp b/src/storm/solver/LinearEquationSolver.cpp
index 193f2cc09..44eb42db0 100644
--- a/src/storm/solver/LinearEquationSolver.cpp
+++ b/src/storm/solver/LinearEquationSolver.cpp
@@ -31,89 +31,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void LinearEquationSolver<ValueType>::repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n) const {
-            if (!cachedRowVector) {
-                cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
-            }
-            
-            // We enable caching for this. But remember how the old setting was
-            bool cachingWasEnabled = isCachingEnabled();
-            setCachingEnabled(true);
-            
-            // Set up some temporary variables so that we can just swap pointers instead of copying the result after
-            // each iteration.
-            std::vector<ValueType>* currentX = &x;
-            std::vector<ValueType>* nextX = cachedRowVector.get();
-            
-            // Now perform matrix-vector multiplication as long as we meet the bound.
-            this->startMeasureProgress();
-            for (uint_fast64_t i = 0; i < n; ++i) {
-                this->multiply(*currentX, b, *nextX);
-                std::swap(nextX, currentX);
-
-                // Potentially show progress.
-                this->showProgressIterative(i, n);
-            }
-            
-            // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x,
-            // because the output is supposed to be stored in the input vector x.
-            if (currentX == cachedRowVector.get()) {
-                std::swap(x, *currentX);
-            }
-            
-            // restore the old caching setting
-            setCachingEnabled(cachingWasEnabled);
-            
-            if (!isCachingEnabled()) {
-                clearCache();
-            }
-        }
-        
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::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) const {
-            if (!cachedRowVector) {
-                cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
-            }
-            
-            // We enable caching for this. But remember how the old setting was
-            bool cachingWasEnabled = isCachingEnabled();
-            setCachingEnabled(true);
-
-            this->multiply(x, b, *cachedRowVector);
-            vectorHelper.reduceVector(dir, *cachedRowVector, result, rowGroupIndices, choices);
-            
-            // restore the old caching setting
-            setCachingEnabled(cachingWasEnabled);
-            
-            if (!isCachingEnabled()) {
-                clearCache();
-            }
-        }
-        
-#ifdef STORM_HAVE_CARL
-        template<>
-        void LinearEquationSolver<storm::RationalFunction>::multiplyAndReduce(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const* b, std::vector<storm::RationalFunction>& result, std::vector<uint_fast64_t>* choices ) const {
-            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Reducing rational function vector is not supported.");
-        }
-#endif
-        
-        template<typename ValueType>
-        bool LinearEquationSolver<ValueType>::supportsGaussSeidelMultiplication() const {
-            return false;
-        }
-        
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::multiplyGaussSeidel(std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
-            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This solver does not support the function 'multiplyGaussSeidel'.");
-        }
-        
-        template<typename ValueType>
-        void LinearEquationSolver<ValueType>::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) const {
-            STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This solver does not support the function 'multiplyAndReduceGaussSeidel'.");
-        }
-        
-        template<typename ValueType>
-        LinearEquationSolverRequirements LinearEquationSolver<ValueType>::getRequirements(Environment const& env, LinearEquationSolverTask const& task) const {
+        LinearEquationSolverRequirements LinearEquationSolver<ValueType>::getRequirements(Environment const& env) const {
             return LinearEquationSolverRequirements();
         }
         
@@ -137,15 +55,15 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::unique_ptr<LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(Environment const& env, storm::storage::SparseMatrix<ValueType> const& matrix, LinearEquationSolverTask const& task) const {
-            std::unique_ptr<LinearEquationSolver<ValueType>> solver = this->create(env, task);
+        std::unique_ptr<LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(Environment const& env, storm::storage::SparseMatrix<ValueType> const& matrix) const {
+            std::unique_ptr<LinearEquationSolver<ValueType>> solver = this->create(env);
             solver->setMatrix(matrix);
             return solver;
         }
         
         template<typename ValueType>
-        std::unique_ptr<LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(Environment const& env, storm::storage::SparseMatrix<ValueType>&& matrix, LinearEquationSolverTask const& task) const {
-            std::unique_ptr<LinearEquationSolver<ValueType>> solver = this->create(env, task);
+        std::unique_ptr<LinearEquationSolver<ValueType>> LinearEquationSolverFactory<ValueType>::create(Environment const& env, storm::storage::SparseMatrix<ValueType>&& matrix) const {
+            std::unique_ptr<LinearEquationSolver<ValueType>> solver = this->create(env);
             solver->setMatrix(std::move(matrix));
             return solver;
         }
@@ -156,8 +74,8 @@ namespace storm {
         }
         
         template<typename ValueType>
-        LinearEquationSolverRequirements LinearEquationSolverFactory<ValueType>::getRequirements(Environment const& env, LinearEquationSolverTask const& task) const {
-            return this->create(env)->getRequirements(env, task);
+        LinearEquationSolverRequirements LinearEquationSolverFactory<ValueType>::getRequirements(Environment const& env) const {
+            return this->create(env)->getRequirements(env);
         }
         
         template<typename ValueType>
@@ -165,9 +83,8 @@ namespace storm {
             // Intentionally left empty.
         }
         
-        
         template<>
-        std::unique_ptr<LinearEquationSolver<storm::RationalNumber>> GeneralLinearEquationSolverFactory<storm::RationalNumber>::create(Environment const& env, LinearEquationSolverTask const& task) const {
+        std::unique_ptr<LinearEquationSolver<storm::RationalNumber>> GeneralLinearEquationSolverFactory<storm::RationalNumber>::create(Environment const& env) const {
             EquationSolverType type = env.solver().getLinearEquationSolverType();
             
              // Adjust the solver type if it is not supported by this value type
@@ -188,7 +105,7 @@ namespace storm {
         }
         
         template<>
-        std::unique_ptr<LinearEquationSolver<storm::RationalFunction>> GeneralLinearEquationSolverFactory<storm::RationalFunction>::create(Environment const& env, LinearEquationSolverTask const& task) const {
+        std::unique_ptr<LinearEquationSolver<storm::RationalFunction>> GeneralLinearEquationSolverFactory<storm::RationalFunction>::create(Environment const& env) const {
             EquationSolverType type = env.solver().getLinearEquationSolverType();
             
              // Adjust the solver type if it is not supported by this value type
@@ -208,11 +125,11 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::unique_ptr<LinearEquationSolver<ValueType>> GeneralLinearEquationSolverFactory<ValueType>::create(Environment const& env, LinearEquationSolverTask const& task) const {
+        std::unique_ptr<LinearEquationSolver<ValueType>> GeneralLinearEquationSolverFactory<ValueType>::create(Environment const& env) const {
             EquationSolverType type = env.solver().getLinearEquationSolverType();
             
             // Adjust the solver type if none was specified and we want sound computations
-            if (env.solver().isForceSoundness() && task != LinearEquationSolverTask::Multiply && type != EquationSolverType::Native && type != EquationSolverType::Eigen && type != EquationSolverType::Elimination && type != EquationSolverType::Topological) {
+            if (env.solver().isForceSoundness() && type != EquationSolverType::Native && type != EquationSolverType::Eigen && type != EquationSolverType::Elimination && type != EquationSolverType::Topological) {
                 if (env.solver().isLinearEquationSolverTypeSetFromDefaultValue()) {
                     type = EquationSolverType::Native;
                     STORM_LOG_INFO("Selecting '" + toString(type) + "' as the linear equation solver to guarantee sound results. If you want to override this, please explicitly specify a different solver.");
diff --git a/src/storm/solver/LinearEquationSolver.h b/src/storm/solver/LinearEquationSolver.h
index cc2f73a7a..c53c23393 100644
--- a/src/storm/solver/LinearEquationSolver.h
+++ b/src/storm/solver/LinearEquationSolver.h
@@ -22,8 +22,7 @@ namespace storm {
     namespace solver {
         
         /*!
-         * An interface that represents an abstract linear equation solver. In addition to solving a system of linear
-         * equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
+         * An interface that represents an abstract linear equation solver.
          */
         template<class ValueType>
         class LinearEquationSolver : public AbstractEquationSolver<ValueType> {
@@ -50,86 +49,6 @@ namespace storm {
              */
             bool solveEquations(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
 
-            /*!
-             * Performs on 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.
-             */
-            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const = 0;
-            
-            /*!
-             * Performs on 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 choices If given, the choices made in the reduction process are written to this vector.
-             */
-            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;
-            
-            /*!
-             * Retrieves whether this solver offers the gauss-seidel style multiplications.
-             */
-            virtual bool supportsGaussSeidelMultiplication() const;
-            
-            /*!
-             * Performs on matrix-vector multiplication x' = A*x + b. It does so in a gauss-seidel style, i.e. reusing
-             * the new x' components in the further multiplication.
-             *
-             * @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.
-             */
-            virtual void multiplyGaussSeidel(std::vector<ValueType>& x, std::vector<ValueType> const* b) const;
-            
-            /*!
-             * Performs on 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. It does so in a gauss-seidel
-             * style, i.e. reusing the new x' components in the further multiplication.
-             *
-             * @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 choices If given, the choices made in the reduction process are written to this vector.
-             */
-            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;
-            
-            /*!
-             * 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
-             * @param offset A value that is added to the matrix-vector multiplication result
-             */
-            virtual ValueType multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x) 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(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n) const;
-            
             /*!
              * Retrieves the format in which this solver expects to solve equations. If the solver expects the equation
              * system format, it solves Ax = b. If it it expects a fixed point format, it solves Ax + b = x.
@@ -140,7 +59,7 @@ namespace storm {
              * Retrieves the requirements of the solver under the current settings. Note that these requirements only
              * apply to solving linear equations and not to the matrix vector multiplications.
              */
-            virtual LinearEquationSolverRequirements getRequirements(Environment const& env, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const;
+            virtual LinearEquationSolverRequirements getRequirements(Environment const& env) const;
             
             /*!
              * Sets whether some of the generated data during solver calls should be cached.
@@ -195,7 +114,7 @@ namespace storm {
              * @param matrix The matrix that defines the equation system.
              * @return A pointer to the newly created solver.
              */
-            std::unique_ptr<LinearEquationSolver<ValueType>> create(Environment const& env, storm::storage::SparseMatrix<ValueType> const& matrix, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const;
+            std::unique_ptr<LinearEquationSolver<ValueType>> create(Environment const& env, storm::storage::SparseMatrix<ValueType> const& matrix) const;
 
             /*!
              * Creates a new linear equation solver instance with the given matrix. The caller gives up posession of the
@@ -204,12 +123,12 @@ namespace storm {
              * @param matrix The matrix that defines the equation system.
              * @return A pointer to the newly created solver.
              */
-            std::unique_ptr<LinearEquationSolver<ValueType>> create(Environment const& env, storm::storage::SparseMatrix<ValueType>&& matrix, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const;
+            std::unique_ptr<LinearEquationSolver<ValueType>> create(Environment const& env, storm::storage::SparseMatrix<ValueType>&& matrix) const;
 
             /*!
              * Creates an equation solver with the current settings, but without a matrix.
              */
-            virtual std::unique_ptr<LinearEquationSolver<ValueType>> create(Environment const& env, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const = 0;
+            virtual std::unique_ptr<LinearEquationSolver<ValueType>> create(Environment const& env) const = 0;
 
             /*!
              * Creates a copy of this factory.
@@ -225,7 +144,7 @@ namespace storm {
              * Retrieves the requirements of the solver if it was created with the current settings. Note that these
              * requirements only apply to solving linear equations and not to the matrix vector multiplications.
              */
-            LinearEquationSolverRequirements getRequirements(Environment const& env, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const;
+            LinearEquationSolverRequirements getRequirements(Environment const& env) const;
         };
 
         template<typename ValueType>
@@ -235,7 +154,7 @@ namespace storm {
             
             using LinearEquationSolverFactory<ValueType>::create;
 
-            virtual std::unique_ptr<LinearEquationSolver<ValueType>> create(Environment const& env, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const override;
+            virtual std::unique_ptr<LinearEquationSolver<ValueType>> create(Environment const& env) const override;
 
             virtual std::unique_ptr<LinearEquationSolverFactory<ValueType>> clone() const override;
         };
diff --git a/src/storm/solver/TopologicalLinearEquationSolver.cpp b/src/storm/solver/TopologicalLinearEquationSolver.cpp
index 6778970c9..e3e733407 100644
--- a/src/storm/solver/TopologicalLinearEquationSolver.cpp
+++ b/src/storm/solver/TopologicalLinearEquationSolver.cpp
@@ -138,14 +138,13 @@ namespace storm {
             if (hasDiagonalEntry) {
                 xi /= denominator;
             }
-            //std::cout << "Solved trivial scc " << sccState << " with result " << globalX[sccState] << std::endl;
             return true;
         }
         
         template<typename ValueType>
         bool TopologicalLinearEquationSolver<ValueType>::solveFullyConnectedEquationSystem(storm::Environment const& sccSolverEnvironment, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             if (!this->sccSolver) {
-                this->sccSolver = GeneralLinearEquationSolverFactory<ValueType>().create(sccSolverEnvironment, LinearEquationSolverTask::SolveEquations);
+                this->sccSolver = GeneralLinearEquationSolverFactory<ValueType>().create(sccSolverEnvironment);
                 this->sccSolver->setCachingEnabled(true);
                 this->sccSolver->setBoundsFromOtherSolver(*this);
                 if (this->sccSolver->getEquationProblemFormat(sccSolverEnvironment) == LinearEquationSolverProblemFormat::EquationSystem) {
@@ -165,7 +164,7 @@ namespace storm {
             
             // Set up the SCC solver
             if (!this->sccSolver) {
-                this->sccSolver = GeneralLinearEquationSolverFactory<ValueType>().create(sccSolverEnvironment, LinearEquationSolverTask::SolveEquations);
+                this->sccSolver = GeneralLinearEquationSolverFactory<ValueType>().create(sccSolverEnvironment);
                 this->sccSolver->setCachingEnabled(true);
             }
             
@@ -215,75 +214,15 @@ namespace storm {
             return returnvalue;
         }
         
-        
-        template<typename ValueType>
-        void TopologicalLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
-            if (&x != &result) {
-                multiplier.multAdd(*A, x, b, result);
-            } else {
-                // If the two vectors are aliases, we need to create a temporary.
-                if (!this->cachedRowVector) {
-                    this->cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
-                }
-                
-                multiplier.multAdd(*A, x, b, *this->cachedRowVector);
-                result.swap(*this->cachedRowVector);
-                
-                if (!this->isCachingEnabled()) {
-                    clearCache();
-                }
-            }
-        }
-        
-        template<typename ValueType>
-        void TopologicalLinearEquationSolver<ValueType>::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) const {
-            if (&x != &result) {
-                multiplier.multAddReduce(dir, rowGroupIndices, *A, x, b, result, choices);
-            } else {
-                // If the two vectors are aliases, we need to create a temporary.
-                if (!this->cachedRowVector) {
-                    this->cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount());
-                }
-            
-                multiplier.multAddReduce(dir, rowGroupIndices, *A, x, b, *this->cachedRowVector, choices);
-                result.swap(*this->cachedRowVector);
-                
-                if (!this->isCachingEnabled()) {
-                    clearCache();
-                }
-            }
-        }
-        
-        template<typename ValueType>
-        bool TopologicalLinearEquationSolver<ValueType>::supportsGaussSeidelMultiplication() const {
-            return true;
-        }
-        
-        template<typename ValueType>
-        void TopologicalLinearEquationSolver<ValueType>::multiplyGaussSeidel(std::vector<ValueType>& x, std::vector<ValueType> const* b) const {
-            STORM_LOG_ASSERT(this->A->getRowCount() == this->A->getColumnCount(), "This function is only applicable for square matrices.");
-            multiplier.multAddGaussSeidelBackward(*A, x, b);
-        }
-        
-        template<typename ValueType>
-        void TopologicalLinearEquationSolver<ValueType>::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) const {
-            multiplier.multAddReduceGaussSeidelBackward(dir, rowGroupIndices, *A, x, b, choices);
-        }
-        
-        template<typename ValueType>
-        ValueType TopologicalLinearEquationSolver<ValueType>::multiplyRow(uint64_t const& rowIndex, std::vector<ValueType> const& x) const {
-            return multiplier.multiplyRow(*A, rowIndex, x);
-        }
-        
         template<typename ValueType>
         LinearEquationSolverProblemFormat TopologicalLinearEquationSolver<ValueType>::getEquationProblemFormat(Environment const& env) const {
             return LinearEquationSolverProblemFormat::FixedPointSystem;
         }
         
         template<typename ValueType>
-        LinearEquationSolverRequirements TopologicalLinearEquationSolver<ValueType>::getRequirements(Environment const& env, LinearEquationSolverTask const& task) const {
+        LinearEquationSolverRequirements TopologicalLinearEquationSolver<ValueType>::getRequirements(Environment const& env) const {
             // Return the requirements of the underlying solver
-            return GeneralLinearEquationSolverFactory<ValueType>().getRequirements(getEnvironmentForUnderlyingSolver(env), task);
+            return GeneralLinearEquationSolverFactory<ValueType>().getRequirements(getEnvironmentForUnderlyingSolver(env));
         }
         
         template<typename ValueType>
@@ -305,7 +244,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> TopologicalLinearEquationSolverFactory<ValueType>::create(Environment const& env, LinearEquationSolverTask const& task) const {
+        std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> TopologicalLinearEquationSolverFactory<ValueType>::create(Environment const& env) const {
             return std::make_unique<storm::solver::TopologicalLinearEquationSolver<ValueType>>();
         }
         
diff --git a/src/storm/solver/TopologicalLinearEquationSolver.h b/src/storm/solver/TopologicalLinearEquationSolver.h
index fcf630706..5a36ec9af 100644
--- a/src/storm/solver/TopologicalLinearEquationSolver.h
+++ b/src/storm/solver/TopologicalLinearEquationSolver.h
@@ -22,15 +22,8 @@ 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, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const override;
+            virtual LinearEquationSolverRequirements getRequirements(Environment const& env) const override;
 
             virtual void clearCache() const override;
 
@@ -77,7 +70,7 @@ namespace storm {
         public:
             using LinearEquationSolverFactory<ValueType>::create;
             
-            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(Environment const& env, LinearEquationSolverTask const& task = LinearEquationSolverTask::Unspecified) const override;
+            virtual std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> create(Environment const& env) const override;
             
             virtual std::unique_ptr<LinearEquationSolverFactory<ValueType>> clone() const override;