diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
index 4ccf61e64..dab1449de 100644
--- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
+++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
@@ -627,19 +627,19 @@ namespace storm {
                         result = std::vector<ValueType>(values.size());
                     }
                 }
-                std::vector<ValueType> multiplicationResult(result.size());
                 
                 std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(uniformizedMatrix));
+                solver->allocateAuxMemory(storm::solver::LinearEquationSolverOperation::MultiplyRepeatedly);
                 
                 if (!useMixedPoissonProbabilities && std::get<0>(foxGlynnResult) > 1) {
                     // Perform the matrix-vector multiplications (without adding).
-                    solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1, &multiplicationResult);
+                    solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1);
                 } else if (useMixedPoissonProbabilities) {
                     std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; };
                     
                     // For the iterations below the left truncation point, we need to add and scale the result with the uniformization rate.
                     for (uint_fast64_t index = 1; index < startingIteration; ++index) {
-                        solver->repeatedMultiply(values, nullptr, 1, &multiplicationResult);
+                        solver->repeatedMultiply(values, nullptr, 1);
                         storm::utility::vector::applyPointwise(result, values, result, addAndScale);
                     }
                 }
@@ -649,7 +649,7 @@ namespace storm {
                 ValueType weight = 0;
                 std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; };
                 for (uint_fast64_t index = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) {
-                    solver->repeatedMultiply(values, addVector, 1, &multiplicationResult);
+                    solver->repeatedMultiply(values, addVector, 1);
                     
                     weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
                     storm::utility::vector::applyPointwise(result, values, result, addAndScale);
@@ -658,7 +658,6 @@ namespace storm {
                 return result;
             }
             
-            
             template <typename ValueType>
             storm::storage::SparseMatrix<ValueType> SparseCtmcCslHelper::computeProbabilityMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) {
                 // Turn the rates into probabilities by scaling each row with the exit rate of the state.
diff --git a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index e56b21d66..7817dc644 100644
--- a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -86,6 +86,7 @@ namespace storm {
                 }
                 
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(aProbabilistic);
+                solver->allocateAuxMemory(storm::solver::MinMaxLinearEquationSolverOperation::SolveEquations);
                 
                 // Perform the actual value iteration
                 // * loop until the step bound has been reached
@@ -94,15 +95,13 @@ namespace storm {
                 //      and 1_G being the characteristic vector for all goal states.
                 // *    perform one timed-step using v_MS := A_MSwG * v_MS + A_MStoPS * v_PS + (A * 1_G)|MS
                 std::vector<ValueType> markovianNonGoalValuesSwap(markovianNonGoalValues);
-                std::vector<ValueType> multiplicationResultScratchMemory(aProbabilistic.getRowCount());
-                std::vector<ValueType> aProbabilisticScratchMemory(probabilisticNonGoalValues.size());
                 for (uint_fast64_t currentStep = 0; currentStep < numberOfSteps; ++currentStep) {
                     // Start by (re-)computing bProbabilistic = bProbabilisticFixed + aProbabilisticToMarkovian * vMarkovian.
                     aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
                     storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
                     
                     // Now perform the inner value iteration for probabilistic states.
-                    solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
+                    solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic);
                     
                     // (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic.
                     aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian);
@@ -115,7 +114,7 @@ namespace storm {
                 // After the loop, perform one more step of the value iteration for PS states.
                 aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
                 storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
-                solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
+                solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic);
             }
 
             template<typename ValueType>
diff --git a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
index 4863c1daf..2f3111a7f 100644
--- a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
@@ -141,7 +141,7 @@ namespace storm {
                     std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
                     
                     std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first));
-                    solver->multiply(dir, x, &explicitRepresentation.second, stepBound);
+                    solver->repeatedMultiply(dir, x, &explicitRepresentation.second, stepBound);
                     
                     // Return a hybrid check result that stores the numerical values explicitly.
                     return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.template toAdd<ValueType>(), maybeStates, odd, x));
@@ -166,7 +166,7 @@ namespace storm {
                 
                 // Perform the matrix-vector multiplication.
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitMatrix));
-                solver->multiply(dir, x, nullptr, stepBound);
+                solver->repeatedMultiply(dir, x, nullptr, stepBound);
                 
                 // Return a hybrid check result that stores the numerical values explicitly.
                 return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero<ValueType>(), model.getReachableStates(), odd, x));
@@ -195,7 +195,7 @@ namespace storm {
                 
                 // Perform the matrix-vector multiplication.
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first));
-                solver->multiply(dir, x, &explicitRepresentation.second, stepBound);
+                solver->repeatedMultiply(dir, x, &explicitRepresentation.second, stepBound);
                 
                 // Return a hybrid check result that stores the numerical values explicitly.
                 return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero<ValueType>(), model.getReachableStates(), odd, x));
diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
index 916d24d53..79a85fdfe 100644
--- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
@@ -121,7 +121,7 @@ namespace storm {
                 
                 // Perform one single matrix-vector multiplication.
                 std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix);
-                solver->repeatedMultiply(result);
+                solver->repeatedMultiply(result, nullptr, 1);
                 return result;
             }
             
diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
index 8c09b96b8..75fbfdfea 100644
--- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
+++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
@@ -49,7 +49,7 @@ namespace storm {
                     std::vector<ValueType> subresult(maybeStates.getNumberOfSetBits());
                     
                     std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(std::move(submatrix));
-                    solver->multiply(dir, subresult, &b, stepBound);
+                    solver->repeatedMultiply(dir, subresult, &b, stepBound);
                     
                     // Set the values of the resulting vector accordingly.
                     storm::utility::vector::setVectorValues(result, maybeStates, subresult);
@@ -67,7 +67,7 @@ namespace storm {
                 storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
                 
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
-                solver->multiply(dir, result);
+                solver->repeatedMultiply(dir, result, nullptr, 1);
                 
                 return result;
             }
@@ -211,7 +211,7 @@ namespace storm {
                 std::vector<ValueType> result(rewardModel.getStateRewardVector());
                 
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
-                solver->multiply(dir, result, nullptr, stepCount);
+                solver->repeatedMultiply(dir, result, nullptr, stepCount);
                 
                 return result;
             }
@@ -235,7 +235,7 @@ namespace storm {
                 }
                 
                 std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
-                solver->multiply(dir, result, &totalRewardVector, stepBound);
+                solver->repeatedMultiply(dir, result, &totalRewardVector, stepBound);
                 
                 return result;
             }
@@ -601,13 +601,13 @@ namespace storm {
                     if (fixedTargetStates.get(state)) {
                         builder.addNextValue(currentRow, newGoalState, conditionProbabilities[state]);
                         if (!storm::utility::isZero(conditionProbabilities[state])) {
-                            builder.addNextValue(currentRow, newFailState, 1 - conditionProbabilities[state]);
+                            builder.addNextValue(currentRow, newFailState, storm::utility::one<ValueType>() - conditionProbabilities[state]);
                         }
                         ++currentRow;
                     } else if (conditionStates.get(state)) {
                         builder.addNextValue(currentRow, newGoalState, targetProbabilities[state]);
                         if (!storm::utility::isZero(targetProbabilities[state])) {
-                            builder.addNextValue(currentRow, newStopState, 1 - targetProbabilities[state]);
+                            builder.addNextValue(currentRow, newStopState, storm::utility::one<ValueType>() - targetProbabilities[state]);
                         }
                         ++currentRow;
                     } else {
diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp
index c00fe97b0..25c16f1ae 100644
--- a/src/solver/EigenLinearEquationSolver.cpp
+++ b/src/solver/EigenLinearEquationSolver.cpp
@@ -109,12 +109,23 @@ namespace storm {
 
         template<typename ValueType>
         EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EigenLinearEquationSolverSettings<ValueType> const& settings) : settings(settings) {
+            this->setMatrix(std::move(A));
+        }
+        
+        template<typename ValueType>
+        void EigenLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
+            eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A);
+        }
+        
+        template<typename ValueType>
+        void EigenLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
+            // Take ownership of the matrix so it is destroyed after we have translated it to Eigen's format.
             storm::storage::SparseMatrix<ValueType> localA(std::move(A));
-            eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(localA);
+            this->setMatrix(localA);
         }
         
         template<typename ValueType>
-        void EigenLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+        void EigenLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             // Map the input vectors to Eigen's format.
             auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             auto eigenB = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b.data(), b.size());
@@ -197,7 +208,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void EigenLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
+        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(Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
 
@@ -234,10 +245,20 @@ namespace storm {
             return settings;
         }
         
+        template<typename ValueType>
+        uint64_t EigenLinearEquationSolver<ValueType>::getMatrixRowCount() const {
+            return eigenA->rows();
+        }
+        
+        template<typename ValueType>
+        uint64_t EigenLinearEquationSolver<ValueType>::getMatrixColumnCount() const {
+            return eigenA->cols();
+        }
+        
         // Specialization form storm::RationalNumber
         
         template<>
-        void EigenLinearEquationSolver<storm::RationalNumber>::solveEquations(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>* multiplyResult) const {
+        void EigenLinearEquationSolver<storm::RationalNumber>::solveEquations(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b) const {
             // Map the input vectors to Eigen's format.
             auto eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             auto eigenB = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b.data(), b.size());
@@ -250,7 +271,7 @@ namespace storm {
         // Specialization form storm::RationalFunction
         
         template<>
-        void EigenLinearEquationSolver<storm::RationalFunction>::solveEquations(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b, std::vector<storm::RationalFunction>* multiplyResult) const {
+        void EigenLinearEquationSolver<storm::RationalFunction>::solveEquations(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b) const {
             // Map the input vectors to Eigen's format.
             auto eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(x.data(), x.size());
             auto eigenB = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b.data(), b.size());
diff --git a/src/solver/EigenLinearEquationSolver.h b/src/solver/EigenLinearEquationSolver.h
index 43e34427e..04c6c9718 100644
--- a/src/solver/EigenLinearEquationSolver.h
+++ b/src/solver/EigenLinearEquationSolver.h
@@ -61,13 +61,19 @@ namespace storm {
             EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EigenLinearEquationSolverSettings<ValueType> const& settings = EigenLinearEquationSolverSettings<ValueType>());
             EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EigenLinearEquationSolverSettings<ValueType> const& settings = EigenLinearEquationSolverSettings<ValueType>());
             
-            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
+            
+            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
+            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
 
             EigenLinearEquationSolverSettings<ValueType>& getSettings();
             EigenLinearEquationSolverSettings<ValueType> const& getSettings() const;
-            
+                        
         private:
+            virtual uint64_t getMatrixRowCount() const override;
+            virtual uint64_t getMatrixColumnCount() const override;
+            
             // The (eigen) matrix associated with this equation solver.
             std::unique_ptr<Eigen::SparseMatrix<ValueType>> eigenA;
 
diff --git a/src/solver/EliminationLinearEquationSolver.cpp b/src/solver/EliminationLinearEquationSolver.cpp
index 75340d6ba..982d2a54b 100644
--- a/src/solver/EliminationLinearEquationSolver.cpp
+++ b/src/solver/EliminationLinearEquationSolver.cpp
@@ -34,19 +34,29 @@ namespace storm {
         }
         
         template<typename ValueType>
-        EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), settings(settings) {
-            // Intentionally left empty.
+        EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings) {
+            this->setMatrix(A);
         }
         
         template<typename ValueType>
-        EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), settings(settings) {
-            // Intentionally left empty.
+        EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings) {
+            this->setMatrix(std::move(A));
         }
         
         template<typename ValueType>
-        void EliminationLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
-            STORM_LOG_WARN_COND(multiplyResult == nullptr, "Providing scratch memory will not improve the performance of this solver.");
-            
+        void EliminationLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
+            this->A = &A;
+            localA.reset();
+        }
+        
+        template<typename ValueType>
+        void EliminationLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
+            localA = std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A));
+            this->A = localA.get();
+        }
+        
+        template<typename ValueType>
+        void EliminationLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             // FIXME: This solver will not work for all input systems. More concretely, the current implementation will
             // not work for systems that have a 0 on the diagonal. This is not a restriction of this technique in general
             // but arbitrary matrices require pivoting, which is not currently implemented.
@@ -59,7 +69,7 @@ namespace storm {
             if (localA) {
                 localA->convertToEquationSystem();
             } else {
-                locallyConvertedMatrix = A;
+                locallyConvertedMatrix = *A;
                 locallyConvertedMatrix.convertToEquationSystem();
             }
             
@@ -101,16 +111,16 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void EliminationLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
+        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);
+                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);
+                A->multiplyWithVector(x, tmp);
                 if (b != nullptr) {
                     storm::utility::vector::addVectors(tmp, *b, result);
                 }
@@ -127,6 +137,16 @@ namespace storm {
             return settings;
         }
         
+        template<typename ValueType>
+        uint64_t EliminationLinearEquationSolver<ValueType>::getMatrixRowCount() const {
+            return this->A->getRowCount();
+        }
+        
+        template<typename ValueType>
+        uint64_t EliminationLinearEquationSolver<ValueType>::getMatrixColumnCount() const {
+            return this->A->getColumnCount();
+        }
+        
         template<typename ValueType>
         std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EliminationLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
             return std::make_unique<storm::solver::EliminationLinearEquationSolver<ValueType>>(matrix, settings);
diff --git a/src/solver/EliminationLinearEquationSolver.h b/src/solver/EliminationLinearEquationSolver.h
index 15a1e5c25..5963c5763 100644
--- a/src/solver/EliminationLinearEquationSolver.h
+++ b/src/solver/EliminationLinearEquationSolver.h
@@ -29,8 +29,11 @@ namespace storm {
             EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EliminationLinearEquationSolverSettings<ValueType> const& settings = EliminationLinearEquationSolverSettings<ValueType>());
             EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings = EliminationLinearEquationSolverSettings<ValueType>());
             
-            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
+            
+            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
+            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
 
             EliminationLinearEquationSolverSettings<ValueType>& getSettings();
             EliminationLinearEquationSolverSettings<ValueType> const& getSettings() const;
@@ -38,13 +41,16 @@ namespace storm {
         private:
             void initializeSettings();
             
+            virtual uint64_t getMatrixRowCount() const override;
+            virtual uint64_t getMatrixColumnCount() const override;
+            
             // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
             // when the solver is destructed.
             std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
 
-            // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
-            // the reference refers to localA.
-            storm::storage::SparseMatrix<ValueType> const& A;
+            // A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix
+            // the pointer refers to localA.
+            storm::storage::SparseMatrix<ValueType> const* A;
             
             // The settings used by the solver.
             EliminationLinearEquationSolverSettings<ValueType> settings;
diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp
index d29cda389..76d38f726 100644
--- a/src/solver/GmmxxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxLinearEquationSolver.cpp
@@ -111,17 +111,31 @@ namespace storm {
         }
         
         template<typename ValueType>
-        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), settings(settings) {
-            // Intentionally left empty.
+        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiMemory(nullptr) {
+            this->setMatrix(A);
         }
 
         template<typename ValueType>
-        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(*localA)), settings(settings) {
-            // Intentionally left empty.
+        GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiMemory(nullptr) {
+            this->setMatrix(std::move(A));
         }
         
         template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+        void GmmxxLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
+            localA.reset();
+            this->A = &A;
+            gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
+            localA = std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A));
+            this->A = localA.get();
+            gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(*localA);
+        }
+        
+        template<typename ValueType>
+        void GmmxxLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             auto method = this->getSettings().getSolutionMethod();
             auto preconditioner = this->getSettings().getPreconditioner();
             STORM_LOG_INFO("Using method '" << method << "' with preconditioner '" << preconditioner << "' (max. " << this->getSettings().getMaximalNumberOfIterations() << " iterations).");
@@ -166,7 +180,7 @@ namespace storm {
                     STORM_LOG_WARN("Iterative solver did not converge.");
                 }
             } else if (method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
-                uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b, multiplyResult);
+                uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(*A, x, b);
                 
                 // Check if the solver converged and issue a warning otherwise.
                 if (iterations < this->getSettings().getMaximalNumberOfIterations()) {
@@ -178,7 +192,7 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void GmmxxLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
+        void GmmxxLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
             if (b) {
                 gmm::mult_add(*gmmxxMatrix, x, *b, result);
             } else {
@@ -187,25 +201,20 @@ namespace storm {
         }
         
         template<typename ValueType>
-        uint_fast64_t GmmxxLinearEquationSolver<ValueType>::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+        uint_fast64_t GmmxxLinearEquationSolver<ValueType>::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+            bool allocatedAuxMemory = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations);
+            if (allocatedAuxMemory) {
+                this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
+            }
+            
             // Get a Jacobi decomposition of the matrix A.
             std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
             
             // Convert the LU matrix to gmm++'s format.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.first));
         
-            // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
-            bool multiplyResultProvided = true;
-            std::vector<ValueType>* nextX = multiplyResult;
-            if (nextX == nullptr) {
-                nextX = new std::vector<ValueType>(x.size());
-                multiplyResultProvided = false;
-            }
-            std::vector<ValueType> const* copyX = nextX;
             std::vector<ValueType>* currentX = &x;
-            
-            // Target vector for precision calculation.
-            std::vector<ValueType> tmpX(x.size());
+            std::vector<ValueType>* nextX = auxiliaryJacobiMemory.get();
             
             // Set up additional environment variables.
             uint_fast64_t iterationCount = 0;
@@ -213,9 +222,8 @@ namespace storm {
             
             while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
                 // Compute D^-1 * (b - LU * x) and store result in nextX.
-                gmm::mult(*gmmLU, *currentX, tmpX);
-                gmm::add(b, gmm::scaled(tmpX, -storm::utility::one<ValueType>()), tmpX);
-                storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
+                gmm::mult_add(*gmmLU, gmm::scaled(*currentX, -storm::utility::one<ValueType>()), b, *nextX);
+                storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX);
                 
                 // Now check if the process already converged within our precision.
                 converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion());
@@ -229,13 +237,13 @@ namespace storm {
             
             // If the last iteration did not write to the original x we have to swap the contents, because the
             // output has to be written to the input parameter x.
-            if (currentX == copyX) {
+            if (currentX == auxiliaryJacobiMemory.get()) {
                 std::swap(x, *currentX);
             }
             
-            // If the vector for the temporary multiplication result was not provided, we need to delete it.
-            if (!multiplyResultProvided) {
-                delete copyX;
+            // If we allocated auxiliary memory, we need to dispose of it now.
+            if (allocatedAuxMemory) {
+                this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
             }
             
             return iterationCount;
@@ -251,6 +259,67 @@ namespace storm {
             return settings;
         }
         
+        template<typename ValueType>
+        bool GmmxxLinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::SolveEquations) {
+                if (this->getSettings().getSolutionMethod() == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
+                    if (!auxiliaryJacobiMemory) {
+                        auxiliaryJacobiMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount());
+                        result = true;
+                    }
+                }
+            }
+            result |= LinearEquationSolver<ValueType>::allocateAuxMemory(operation);
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool GmmxxLinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::SolveEquations) {
+                if (auxiliaryJacobiMemory) {
+                    result = true;
+                    auxiliaryJacobiMemory.reset();
+                }
+            }
+            result |= LinearEquationSolver<ValueType>::deallocateAuxMemory(operation);
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool GmmxxLinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::SolveEquations) {
+                if (auxiliaryJacobiMemory) {
+                    result = auxiliaryJacobiMemory->size() != this->getMatrixColumnCount();
+                    auxiliaryJacobiMemory->resize(this->getMatrixRowCount());
+                }
+            }
+            result |= LinearEquationSolver<ValueType>::reallocateAuxMemory(operation);
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool GmmxxLinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::SolveEquations) {
+                result |= static_cast<bool>(auxiliaryJacobiMemory);
+            }
+            result |= LinearEquationSolver<ValueType>::hasAuxMemory(operation);
+            return result;
+        }
+        
+        template<typename ValueType>
+        uint64_t GmmxxLinearEquationSolver<ValueType>::getMatrixRowCount() const {
+            return this->A->getRowCount();
+        }
+        
+        template<typename ValueType>
+        uint64_t GmmxxLinearEquationSolver<ValueType>::getMatrixColumnCount() const {
+            return this->A->getColumnCount();
+        }
+        
         template<typename ValueType>
         std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
             return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(matrix, settings);
diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h
index 4425ec940..e7ab4c828 100644
--- a/src/solver/GmmxxLinearEquationSolver.h
+++ b/src/solver/GmmxxLinearEquationSolver.h
@@ -89,12 +89,20 @@ namespace storm {
             GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings = GmmxxLinearEquationSolverSettings<ValueType>());
             GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings = GmmxxLinearEquationSolverSettings<ValueType>());
             
-            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
+            
+            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
+            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
 
             GmmxxLinearEquationSolverSettings<ValueType>& getSettings();
             GmmxxLinearEquationSolverSettings<ValueType> const& getSettings() const;
             
+            virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const override;
+            virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const override;
+            virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const override;
+            virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const override;
+            
         private:
             /*!
              * Solves the linear equation system A*x = b given by the parameters using the Jacobi method.
@@ -108,21 +116,27 @@ namespace storm {
              * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
              * vector must be equal to the number of rows of A.
              */
-            uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const;
+            uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            
+            virtual uint64_t getMatrixRowCount() const override;
+            virtual uint64_t getMatrixColumnCount() const override;
 
             // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
             // when the solver is destructed.
             std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
 
-            // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
-            // the reference refers to localA.
-            storm::storage::SparseMatrix<ValueType> const& A;
+            // A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix
+            // the pointer refers to localA.
+            storm::storage::SparseMatrix<ValueType> const* A;
             
             // The (gmm++) matrix associated with this equation solver.
             std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix;
             
             // The settings used by the solver.
             GmmxxLinearEquationSolverSettings<ValueType> settings;
+            
+            // Auxiliary storage for the Jacobi method.
+            mutable std::unique_ptr<std::vector<ValueType>> auxiliaryJacobiMemory;
         };
         
         template<typename ValueType>
diff --git a/src/solver/LinearEquationSolver.cpp b/src/solver/LinearEquationSolver.cpp
index 84f5f84e5..81dc4124b 100644
--- a/src/solver/LinearEquationSolver.cpp
+++ b/src/solver/LinearEquationSolver.cpp
@@ -14,36 +14,78 @@ namespace storm {
     namespace solver {
         
         template<typename ValueType>
-        void LinearEquationSolver<ValueType>::repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+        LinearEquationSolver<ValueType>::LinearEquationSolver() : auxiliaryRepeatedMultiplyMemory(nullptr) {
+            // Intentionally left empty.
+        }
+        
+        template<typename ValueType>
+        void LinearEquationSolver<ValueType>::repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n) const {
+            bool allocatedAuxMemory = !this->hasAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly);
+            if (allocatedAuxMemory) {
+                this->allocateAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly);
+            }
             
             // 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;
-            
-            bool multiplyResultProvided = true;
-            std::vector<ValueType>* nextX = multiplyResult;
-            if (nextX == nullptr) {
-                nextX = new std::vector<ValueType>(x.size());
-                multiplyResultProvided = false;
-            }
-            std::vector<ValueType> const* copyX = nextX;
+            std::vector<ValueType>* nextX = auxiliaryRepeatedMultiplyMemory.get();
             
             // Now perform matrix-vector multiplication as long as we meet the bound.
             for (uint_fast64_t i = 0; i < n; ++i) {
-                this->multiply(*currentX, *nextX, b);
+                this->multiply(*currentX, b, *nextX);
                 std::swap(nextX, currentX);
             }
             
             // 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 == copyX) {
+            if (currentX == auxiliaryRepeatedMultiplyMemory.get()) {
                 std::swap(x, *currentX);
             }
-            
-            // If the vector for the temporary multiplication result was not provided, we need to delete it.
-            if (!multiplyResultProvided) {
-                delete copyX;
+
+            // If we allocated auxiliary memory, we need to dispose of it now.
+            if (allocatedAuxMemory) {
+                this->deallocateAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly);
+            }
+        }
+        
+        template<typename ValueType>
+        bool LinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const {
+            if (!auxiliaryRepeatedMultiplyMemory) {
+                auxiliaryRepeatedMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixColumnCount());
+                return true;
+            }
+            return false;
+        }
+        
+        template<typename ValueType>
+        bool LinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const {
+            if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) {
+                if (auxiliaryRepeatedMultiplyMemory) {
+                    auxiliaryRepeatedMultiplyMemory.reset();
+                    return true;
+                }
+            }
+            return false;
+        }
+        
+        template<typename ValueType>
+        bool LinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) {
+                if (auxiliaryRepeatedMultiplyMemory) {
+                    result = auxiliaryRepeatedMultiplyMemory->size() != this->getMatrixColumnCount();
+                    auxiliaryRepeatedMultiplyMemory->resize(this->getMatrixColumnCount());
+                }
+            }
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool LinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const {
+            if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) {
+                return static_cast<bool>(auxiliaryRepeatedMultiplyMemory);
             }
+            return false;
         }
         
         template<typename ValueType>
diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h
index 9632ddb2b..aa243d593 100644
--- a/src/solver/LinearEquationSolver.h
+++ b/src/solver/LinearEquationSolver.h
@@ -11,6 +11,10 @@
 namespace storm {
     namespace solver {
         
+        enum class LinearEquationSolverOperation {
+            SolveEquations, MultiplyRepeatedly
+        };
+        
         /*!
          * 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.
@@ -18,10 +22,15 @@ namespace storm {
         template<class ValueType>
         class LinearEquationSolver : public AbstractEquationSolver<ValueType> {
         public:
+            LinearEquationSolver();
+            
             virtual ~LinearEquationSolver() {
                 // Intentionally left empty.
             }
             
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) = 0;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) = 0;
+            
             /*!
              * Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution.
              * The solution of the set of linear equations will be written to the vector x. Note that the matrix A has
@@ -29,22 +38,20 @@ namespace storm {
              *
              * @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A.
              * @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A.
-             * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
-             * vector must be equal to the number of rows of A.
              */
-            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
+            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const = 0;
             
             /*!
              * 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 result The target vector into which to write the multiplication result. Its length must be equal
-             * to the number of rows 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>& result, std::vector<ValueType> const* b = nullptr) const = 0;
+            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const = 0;
             
             /*!
              * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After
@@ -55,10 +62,57 @@ namespace storm {
              * 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 multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
-             * vector 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;
+            
+            // Methods related to allocating/freeing auxiliary storage.
+            
+            /*!
+             * Allocates auxiliary memory that can be used to perform the provided operation. Repeated calls to the
+             * corresponding function can then be run without allocating/deallocating this memory repeatedly.
+             * Note: Since the allocated memory is fit to the currently selected options of the solver, they must not
+             * be changed any more after allocating the auxiliary memory until it is deallocated again.
+             *
+             * @return True iff auxiliary memory was allocated.
+             */
+            virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const;
+            
+            /*!
+             * Destroys previously allocated auxiliary memory for the provided operation.
+             *
+             * @return True iff auxiliary memory was deallocated.
              */
-            void repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const;
+            virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const;
+            
+            /*!
+             * If the matrix dimensions changed and auxiliary memory was allocated, this function needs to be called to
+             * update the auxiliary memory.
+             *
+             * @return True iff the auxiliary memory was reallocated.
+             */
+            virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const;
+            
+            /*!
+             * Checks whether the solver has allocated auxiliary memory for the provided operation.
+             *
+             * @return True iff auxiliary memory was previously allocated (and not yet deallocated).
+             */
+            virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const;
+            
+        private:
+            /*!
+             * Retrieves the row count of the matrix associated with this solver.
+             */
+            virtual uint64_t getMatrixRowCount() const = 0;
+            
+            /*!
+             * Retrieves the column count of the matrix associated with this solver.
+             */
+            virtual uint64_t getMatrixColumnCount() const = 0;
+            
+            // Auxiliary memory for repeated matrix-vector multiplication.
+            mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRepeatedMultiplyMemory;
         };
         
         template<typename ValueType>
diff --git a/src/solver/MinMaxLinearEquationSolver.cpp b/src/solver/MinMaxLinearEquationSolver.cpp
index bd522834f..a92f19dd9 100644
--- a/src/solver/MinMaxLinearEquationSolver.cpp
+++ b/src/solver/MinMaxLinearEquationSolver.cpp
@@ -28,15 +28,15 @@ namespace storm {
         }
 
         template<typename ValueType>
-        void MinMaxLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+        void MinMaxLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             STORM_LOG_THROW(isSet(this->direction), storm::exceptions::IllegalFunctionCallException, "Optimization direction not set.");
-            solveEquations(convert(this->direction), x, b, multiplyResult, newX);
+            solveEquations(convert(this->direction), x, b);
         }
 
         template<typename ValueType>
-        void MinMaxLinearEquationSolver<ValueType>::multiply( std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
+        void MinMaxLinearEquationSolver<ValueType>::repeatedMultiply( std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const {
             STORM_LOG_THROW(isSet(this->direction), storm::exceptions::IllegalFunctionCallException, "Optimization direction not set.");
-            return multiply(convert(this->direction), x, b, n, multiplyResult);
+            return repeatedMultiply(convert(this->direction), x, b, n);
         }
 
         template<typename ValueType>
@@ -79,6 +79,21 @@ namespace storm {
             return std::move(scheduler.get());
         }
         
+        template<typename ValueType>
+        bool MinMaxLinearEquationSolver<ValueType>::allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
+            return false;
+        }
+        
+        template<typename ValueType>
+        bool MinMaxLinearEquationSolver<ValueType>::deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
+            return false;
+        }
+        
+        template<typename ValueType>
+        bool MinMaxLinearEquationSolver<ValueType>::hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
+            return false;
+        }
+        
         template<typename ValueType>
         MinMaxLinearEquationSolverFactory<ValueType>::MinMaxLinearEquationSolverFactory(bool trackScheduler) : trackScheduler(trackScheduler) {
             // Intentionally left empty.
diff --git a/src/solver/MinMaxLinearEquationSolver.h b/src/solver/MinMaxLinearEquationSolver.h
index 358ef0be1..0582269fb 100644
--- a/src/solver/MinMaxLinearEquationSolver.h
+++ b/src/solver/MinMaxLinearEquationSolver.h
@@ -20,6 +20,10 @@ namespace storm {
     
     namespace solver {
         
+        enum class MinMaxLinearEquationSolverOperation {
+            SolveEquations, MultiplyRepeatedly
+        };
+        
         /*!
          * A class representing the interface that all min-max linear equation solvers shall implement.
          */
@@ -40,20 +44,15 @@ namespace storm {
              * @param x The solution vector x. The initial values of x represent a guess of the real values to the
              * solver, but may be ignored.
              * @param b The vector to add after matrix-vector multiplication.
-             * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
-             * vector must be equal to the number of rows of A.
-             * @param newX If non-null, this memory is used as a scratch memory. If given, the length of this
-             * vector must be equal to the length of the vector x (and thus to the number of columns of A).
-             * @return The solution vector x of the system of linear equations as the content of the parameter x.
              */
-            virtual void solveEquations(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const = 0;
+            virtual void solveEquations(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b) const = 0;
             
             /*!
              * Behaves the same as the other variant of <code>solveEquations</code>, with the distinction that
              * instead of providing the optimization direction as an argument, the internally set optimization direction
              * is used. Note: this method can only be called after setting the optimization direction.
              */
-            void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const;
+            void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             
             /*!
              * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes
@@ -68,18 +67,16 @@ namespace storm {
              * i.e. after the method returns, this vector will contain the computed values.
              * @param b If not null, this vector is added after each multiplication.
              * @param n Specifies the number of iterations the matrix-vector multiplication is performed.
-             * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
-             * vector must be equal to the number of rows of A.
              * @return The result of the repeated matrix-vector multiplication as the content of the vector x.
              */
-            virtual void multiply(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
+            virtual void repeatedMultiply(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n = 1) const = 0;
             
             /*!
              * Behaves the same as the other variant of <code>multiply</code>, with the
              * distinction that instead of providing the optimization direction as an argument, the internally set
              * optimization direction is used. Note: this method can only be called after setting the optimization direction.
              */
-            virtual void multiply( std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const;
+            virtual void repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType>* b , uint_fast64_t n) const;
             
             /*!
              * Sets an optimization direction to use for calls to methods that do not explicitly provide one.
@@ -119,6 +116,32 @@ namespace storm {
              */
             std::unique_ptr<storm::storage::TotalScheduler> getScheduler();
             
+            // Methods related to allocating/freeing auxiliary storage.
+            
+            /*!
+             * Allocates auxiliary storage that can be used to perform the provided operation. Repeated calls to the
+             * corresponding function can then be run without allocating/deallocating this storage repeatedly.
+             * Note: Since the allocated storage is fit to the currently selected options of the solver, they must not
+             * be changed any more after allocating the auxiliary storage until the storage is deallocated again.
+             *
+             * @return True iff auxiliary storage was allocated.
+             */
+            virtual bool allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const;
+            
+            /*!
+             * Destroys previously allocated auxiliary storage for the provided operation.
+             *
+             * @return True iff auxiliary storage was deallocated.
+             */
+            virtual bool deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const;
+                        
+            /*!
+             * Checks whether the solver has allocated auxiliary storage for the provided operation.
+             *
+             * @return True iff auxiliary storage was previously allocated (and not yet deallocated).
+             */
+            virtual bool hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const;
+            
         protected:
             /// The optimization direction to use for calls to functions that do not provide it explicitly. Can also be unset.
             OptimizationDirectionSetting direction;
diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp
index b24c9dde4..ea39a075c 100644
--- a/src/solver/NativeLinearEquationSolver.cpp
+++ b/src/solver/NativeLinearEquationSolver.cpp
@@ -84,70 +84,62 @@ namespace storm {
         }
         
         template<typename ValueType>
-        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), settings(settings) {
-            // Intentionally left empty.
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingMemory(nullptr) {
+            this->setMatrix(A);
         }
 
         template<typename ValueType>
-        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), settings(settings) {
-            // Intentionally left empty.
+        NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingMemory(nullptr) {
+            this->setMatrix(std::move(A));
         }
         
         template<typename ValueType>
-        void NativeLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
+        void NativeLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
+            localA.reset();
+            this->A = &A;
+        }
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
+            localA = std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A));
+            this->A = localA.get();
+        }
+        
+        template<typename ValueType>
+        void NativeLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+            bool allocatedAuxStorage = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations);
+            if (allocatedAuxStorage) {
+                this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
+            }
+            
             if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::GaussSeidel) {
                 // Define the omega used for SOR.
                 ValueType omega = this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::SOR ? this->getSettings().getOmega() : storm::utility::one<ValueType>();
                 
-                // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
-                bool tmpXProvided = true;
-                std::vector<ValueType>* tmpX = multiplyResult;
-                if (multiplyResult == nullptr) {
-                    tmpX = new std::vector<ValueType>(x);
-                    tmpXProvided = false;
-                } else {
-                    *tmpX = x;
-                }
-                
                 // Set up additional environment variables.
                 uint_fast64_t iterationCount = 0;
                 bool converged = false;
                 
                 while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations()) {
-                    A.performSuccessiveOverRelaxationStep(omega, x, b);
+                    A->performSuccessiveOverRelaxationStep(omega, x, b);
                     
                     // Now check if the process already converged within our precision.
-                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(x, *tmpX, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x));
+                    converged = storm::utility::vector::equalModuloPrecision<ValueType>(*auxiliarySolvingMemory, x, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x));
                     
-                    // If we did not yet converge, we need to copy the contents of x to *tmpX.
+                    // If we did not yet converge, we need to backup the contents of x.
                     if (!converged) {
-                        *tmpX = x;
+                        *auxiliarySolvingMemory = x;
                     }
                     
                     // Increase iteration count so we can abort if convergence is too slow.
                     ++iterationCount;
                 }
-                
-                // If the vector for the temporary multiplication result was not provided, we need to delete it.
-                if (!tmpXProvided) {
-                    delete tmpX;
-                }
             } else {
                 // Get a Jacobi decomposition of the matrix A.
-                std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+                std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A->getJacobiDecomposition();
                 
-                // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
-                bool multiplyResultProvided = true;
-                std::vector<ValueType>* nextX = multiplyResult;
-                if (nextX == nullptr) {
-                    nextX = new std::vector<ValueType>(x.size());
-                    multiplyResultProvided = false;
-                }
-                std::vector<ValueType> const* copyX = nextX;
                 std::vector<ValueType>* currentX = &x;
-                
-                // Target vector for multiplication.
-                std::vector<ValueType> tmpX(x.size());
+                std::vector<ValueType>* nextX = auxiliarySolvingMemory.get();
                 
                 // Set up additional environment variables.
                 uint_fast64_t iterationCount = 0;
@@ -155,15 +147,15 @@ namespace storm {
                 
                 while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
                     // Compute D^-1 * (b - LU * x) and store result in nextX.
-                    jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
-                    storm::utility::vector::subtractVectors(b, tmpX, tmpX);
-                    storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
-                    
-                    // Swap the two pointers as a preparation for the next iteration.
-                    std::swap(nextX, currentX);
+                    jacobiDecomposition.first.multiplyWithVector(*currentX, *nextX);
+                    storm::utility::vector::subtractVectors(b, *nextX, *nextX);
+                    storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX);
                     
                     // Now check if the process already converged within our precision.
                     converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion());
+
+                    // Swap the two pointers as a preparation for the next iteration.
+                    std::swap(nextX, currentX);
                     
                     // Increase iteration count so we can abort if convergence is too slow.
                     ++iterationCount;
@@ -171,28 +163,28 @@ namespace storm {
                                 
                 // If the last iteration did not write to the original x we have to swap the contents, because the
                 // output has to be written to the input parameter x.
-                if (currentX == copyX) {
+                if (currentX == auxiliarySolvingMemory.get()) {
                     std::swap(x, *currentX);
                 }
-                
-                // If the vector for the temporary multiplication result was not provided, we need to delete it.
-                if (!multiplyResultProvided) {
-                    delete copyX;
-                }
+            }
+            
+            // If we allocated auxiliary memory, we need to dispose of it now.
+            if (allocatedAuxStorage) {
+                this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
             }
         }
         
         template<typename ValueType>
-        void NativeLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
+        void NativeLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
             if (&x != &result) {
-                A.multiplyWithVector(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);
+                A->multiplyWithVector(x, tmp);
                 if (b != nullptr) {
                     storm::utility::vector::addVectors(tmp, *b, result);
                 }
@@ -209,6 +201,65 @@ namespace storm {
             return settings;
         }
         
+        template<typename ValueType>
+        bool NativeLinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::SolveEquations) {
+                if (!auxiliarySolvingMemory) {
+                    auxiliarySolvingMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount());
+                    result = true;
+                }
+            }
+            result |= LinearEquationSolver<ValueType>::allocateAuxMemory(operation);
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool NativeLinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::SolveEquations) {
+                if (auxiliarySolvingMemory) {
+                    result = true;
+                    auxiliarySolvingMemory.reset();
+                }
+            }
+            result |= LinearEquationSolver<ValueType>::deallocateAuxMemory(operation);
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool NativeLinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::SolveEquations) {
+                if (auxiliarySolvingMemory) {
+                    result = auxiliarySolvingMemory->size() != this->getMatrixColumnCount();
+                    auxiliarySolvingMemory->resize(this->getMatrixRowCount());
+                }
+            }
+            result |= LinearEquationSolver<ValueType>::reallocateAuxMemory(operation);
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool NativeLinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == LinearEquationSolverOperation::SolveEquations) {
+                result |= static_cast<bool>(auxiliarySolvingMemory);
+            }
+            result |= LinearEquationSolver<ValueType>::hasAuxMemory(operation);
+            return result;
+        }
+        
+        template<typename ValueType>
+        uint64_t NativeLinearEquationSolver<ValueType>::getMatrixRowCount() const {
+            return this->A->getRowCount();
+        }
+        
+        template<typename ValueType>
+        uint64_t NativeLinearEquationSolver<ValueType>::getMatrixColumnCount() const {
+            return this->A->getColumnCount();
+        }
+        
         template<typename ValueType>
         std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
             return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(matrix, settings);
diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h
index 235985df6..436e5cd36 100644
--- a/src/solver/NativeLinearEquationSolver.h
+++ b/src/solver/NativeLinearEquationSolver.h
@@ -46,23 +46,37 @@ namespace storm {
             NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings = NativeLinearEquationSolverSettings<ValueType>());
             NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings = NativeLinearEquationSolverSettings<ValueType>());
             
-            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
-            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
+            virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
+            
+            virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
+            virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
 
             NativeLinearEquationSolverSettings<ValueType>& getSettings();
             NativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
 
+            virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const override;
+            virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const override;
+            virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const override;
+            virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const override;
+
         private:
+            virtual uint64_t getMatrixRowCount() const override;
+            virtual uint64_t getMatrixColumnCount() const override;
+
             // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
             // when the solver is destructed.
             std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
             
-            // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
-            // the reference refers to localA.
-            storm::storage::SparseMatrix<ValueType> const& A;
+            // A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix
+            // the pointer refers to localA.
+            storm::storage::SparseMatrix<ValueType> const* A;
                         
             // The settings used by the solver.
             NativeLinearEquationSolverSettings<ValueType> settings;
+            
+            // Auxiliary memory for the equation solving methods.
+            mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingMemory;
         };
         
         template<typename ValueType>
diff --git a/src/solver/StandardMinMaxLinearEquationSolver.cpp b/src/solver/StandardMinMaxLinearEquationSolver.cpp
index db507ef68..4e0609e80 100644
--- a/src/solver/StandardMinMaxLinearEquationSolver.cpp
+++ b/src/solver/StandardMinMaxLinearEquationSolver.cpp
@@ -84,45 +84,38 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             switch (this->getSettings().getSolutionMethod()) {
-                case StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration: solveEquationsValueIteration(dir, x, b, multiplyResult, newX); break;
-                case StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration: solveEquationsPolicyIteration(dir, x, b, multiplyResult, newX); break;
+                case StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration: solveEquationsValueIteration(dir, x, b); break;
+                case StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration: solveEquationsPolicyIteration(dir, x, b); break;
             }
         }
         
         template<typename ValueType>
-        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
-            
-            // Create scratch memory if none was provided.
-            bool multiplyResultMemoryProvided = multiplyResult != nullptr;
-            if (multiplyResult == nullptr) {
-                multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
-            }
-            
+        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             // Create the initial scheduler.
             std::vector<storm::storage::sparse::state_type> scheduler(this->A.getRowGroupCount());
             
             // Create a vector for storing the right-hand side of the inner equation system.
             std::vector<ValueType> subB(this->A.getRowGroupCount());
-            
-            // Create a vector that the inner equation solver can use as scratch memory.
-            std::vector<ValueType> deterministicMultiplyResult(this->A.getRowGroupCount());
+
+            // Resolve the nondeterminism according to the current scheduler.
+            storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(scheduler, true);
+            submatrix.convertToEquationSystem();
+            storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A.getRowGroupIndices(), b);
+
+            // Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration.
+            auto solver = linearEquationSolverFactory->create(std::move(submatrix));
+            solver->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
             
             Status status = Status::InProgress;
             uint64_t iterations = 0;
             do {
-                // Resolve the nondeterminism according to the current scheduler.
-                storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(scheduler, true);
-                submatrix.convertToEquationSystem();
-                storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A.getRowGroupIndices(), b);
-                
                 // Solve the equation system for the 'DTMC'.
                 // FIXME: we need to remove the 0- and 1- states to make the solution unique.
                 // HOWEVER: if we start with a valid scheduler, then we will never get an illegal one, because staying
                 // within illegal MECs will never strictly improve the value. Is this true?
-                auto solver = linearEquationSolverFactory->create(std::move(submatrix));
-                solver->solveEquations(x, subB, &deterministicMultiplyResult);
+                solver->solveEquations(x, subB);
                 
                 // Go through the multiplication result and see whether we can improve any of the choices.
                 bool schedulerImproved = false;
@@ -151,6 +144,12 @@ namespace storm {
                 // If the scheduler did not improve, we are done.
                 if (!schedulerImproved) {
                     status = Status::Converged;
+                } else {
+                    // Update the scheduler and the solver.
+                    submatrix = this->A.selectRowsFromRowGroups(scheduler, true);
+                    submatrix.convertToEquationSystem();
+                    storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A.getRowGroupIndices(), b);
+                    solver->setMatrix(std::move(submatrix));
                 }
                 
                 // Update environment variables.
@@ -164,10 +163,6 @@ namespace storm {
             if (this->isTrackSchedulerSet()) {
                 this->scheduler = std::make_unique<storm::storage::TotalScheduler>(std::move(scheduler));
             }
-            
-            if (!multiplyResultMemoryProvided) {
-                delete multiplyResult;
-            }
         }
         
         template<typename ValueType>
@@ -186,22 +181,15 @@ namespace storm {
         }
         
         template<typename ValueType>
-        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+        void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
-            
-            // Create scratch memory if none was provided.
-            bool multiplyResultMemoryProvided = multiplyResult != nullptr;
-            if (multiplyResult == nullptr) {
-                multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
-            }
-            std::vector<ValueType>* currentX = &x;
-            bool xMemoryProvided = newX != nullptr;
-            if (newX == nullptr) {
-                newX = new std::vector<ValueType>(x.size());
+            bool allocatedAuxMemory = !this->hasAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations);
+            if (allocatedAuxMemory) {
+                this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations);
             }
             
-            // Keep track of which of the vectors for x is the auxiliary copy.
-            std::vector<ValueType>* copyX = newX;
+            std::vector<ValueType>* currentX = &x;
+            std::vector<ValueType>* newX = auxiliarySolvingVectorMemory.get();
             
             // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations.
             uint64_t iterations = 0;
@@ -209,10 +197,10 @@ namespace storm {
             Status status = Status::InProgress;
             while (status == Status::InProgress) {
                 // Compute x' = A*x + b.
-                solver->multiply(*currentX, *multiplyResult, &b);
+                solver->multiply(*currentX, &b, *auxiliarySolvingMultiplyMemory);
                 
                 // Reduce the vector x' by applying min/max for all non-deterministic choices.
-                storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, *newX, this->A.getRowGroupIndices());
+                storm::utility::vector::reduceVectorMinOrMax(dir, *auxiliarySolvingMultiplyMemory, *newX, this->A.getRowGroupIndices());
                 
                 // Determine whether the method converged.
                 if (storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) {
@@ -229,39 +217,36 @@ namespace storm {
             
             // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
             // is currently stored in currentX, but x is the output vector.
-            if (currentX == copyX) {
+            if (currentX == auxiliarySolvingVectorMemory.get()) {
                 std::swap(x, *currentX);
             }
             
-            // Dispose of allocated scratch memory.
-            if (!xMemoryProvided) {
-                delete copyX;
-            }
-            if (!multiplyResultMemoryProvided) {
-                delete multiplyResult;
+            // If we allocated auxiliary memory, we need to dispose of it now.
+            if (allocatedAuxMemory) {
+                this->deallocateAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations);
             }
         }
         
         template<typename ValueType>
-        void StandardMinMaxLinearEquationSolver<ValueType>::multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
-            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
-            
-            // If scratch memory was not provided, we create it.
-            bool multiplyResultMemoryProvided = multiplyResult != nullptr;
-            if (!multiplyResult) {
-                multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
+        void StandardMinMaxLinearEquationSolver<ValueType>::repeatedMultiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const {
+            bool allocatedAuxMemory = !this->hasAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly);
+            if (allocatedAuxMemory) {
+                this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly);
             }
             
+            std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
+            
             for (uint64_t i = 0; i < n; ++i) {
-                solver->multiply(x, *multiplyResult, b);
+                solver->multiply(x, b, *auxiliaryRepeatedMultiplyMemory);
                 
                 // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
                 // element of the min/max operator stack.
-                storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices());
+                storm::utility::vector::reduceVectorMinOrMax(dir, *auxiliaryRepeatedMultiplyMemory, x, this->A.getRowGroupIndices());
             }
             
-            if (!multiplyResultMemoryProvided) {
-                delete multiplyResult;
+            // If we allocated auxiliary memory, we need to dispose of it now.
+            if (allocatedAuxMemory) {
+                this->deallocateAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly);
             }
         }
         
@@ -298,6 +283,73 @@ namespace storm {
             return settings;
         }
         
+        template<typename ValueType>
+        bool StandardMinMaxLinearEquationSolver<ValueType>::allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) {
+                if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) {
+                    if (!auxiliarySolvingMultiplyMemory) {
+                        result = true;
+                        auxiliarySolvingMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowCount());
+                    }
+                    if (!auxiliarySolvingVectorMemory) {
+                        result = true;
+                        auxiliarySolvingVectorMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowGroupCount());
+                    }
+                } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) {
+                    // Nothing to do in this case.
+                } else {
+                    STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method.");
+                }
+            } else {
+                if (!auxiliaryRepeatedMultiplyMemory) {
+                    result = true;
+                    auxiliaryRepeatedMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowCount());
+                }
+            }
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool StandardMinMaxLinearEquationSolver<ValueType>::deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
+            bool result = false;
+            if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) {
+                if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) {
+                    if (auxiliarySolvingMultiplyMemory) {
+                        result = true;
+                        auxiliarySolvingMultiplyMemory.reset();
+                    }
+                    if (auxiliarySolvingVectorMemory) {
+                        result = true;
+                        auxiliarySolvingVectorMemory.reset();
+                    }
+                } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) {
+                    // Nothing to do in this case.
+                } else {
+                    STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method.");
+                }
+            } else {
+                if (auxiliaryRepeatedMultiplyMemory) {
+                    result = true;
+                    auxiliaryRepeatedMultiplyMemory.reset();
+                }
+            }
+            return result;
+        }
+        
+        template<typename ValueType>
+        bool StandardMinMaxLinearEquationSolver<ValueType>::hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
+            if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) {
+                if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) {
+                    return auxiliarySolvingMultiplyMemory && auxiliarySolvingVectorMemory;
+                } else {
+                    return false;
+                }
+            } else {
+                return static_cast<bool>(auxiliaryRepeatedMultiplyMemory);
+            }
+        }
+        
         template<typename ValueType>
         StandardMinMaxLinearEquationSolverFactory<ValueType>::StandardMinMaxLinearEquationSolverFactory(bool trackScheduler) : MinMaxLinearEquationSolverFactory<ValueType>(trackScheduler), linearEquationSolverFactory(nullptr) {
             // Intentionally left empty.
diff --git a/src/solver/StandardMinMaxLinearEquationSolver.h b/src/solver/StandardMinMaxLinearEquationSolver.h
index 36b2aa458..0c0b37f79 100644
--- a/src/solver/StandardMinMaxLinearEquationSolver.h
+++ b/src/solver/StandardMinMaxLinearEquationSolver.h
@@ -38,15 +38,19 @@ namespace storm {
             StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings<ValueType> const& settings = StandardMinMaxLinearEquationSolverSettings<ValueType>());
             StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings<ValueType> const& settings = StandardMinMaxLinearEquationSolverSettings<ValueType>());
             
-            virtual void solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
-            virtual void multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
+            virtual void solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
+            virtual void repeatedMultiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const override;
 
             StandardMinMaxLinearEquationSolverSettings<ValueType> const& getSettings() const;
             StandardMinMaxLinearEquationSolverSettings<ValueType>& getSettings();
             
+            virtual bool allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const override;
+            virtual bool deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const override;
+            virtual bool hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const override;
+            
         private:
-            void solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const;
-            void solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const;
+            void solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            void solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
 
             bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const;
             
@@ -70,6 +74,13 @@ namespace storm {
             // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
             // the reference refers to localA.
             storm::storage::SparseMatrix<ValueType> const& A;
+            
+            // Auxiliary memory for equation solving.
+            mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingMultiplyMemory;
+            mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingVectorMemory;
+            
+            // Auxiliary memory for repeated matrix-vector multiplication.
+            mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRepeatedMultiplyMemory;
         };
      
         template<typename ValueType>
diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
index 6ef7d231e..909c41c30 100644
--- a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
+++ b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
@@ -33,7 +33,7 @@ namespace storm {
         }
                 
         template<typename ValueType>
-		void TopologicalMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
+		void TopologicalMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
 			
 #ifdef GPU_USE_FLOAT
 #define __FORCE_FLOAT_CALCULATION true
@@ -49,7 +49,7 @@ namespace storm {
 				std::vector<float> new_x = storm::utility::vector::toValueType<float>(x);
 				std::vector<float> const new_b = storm::utility::vector::toValueType<float>(b);
 
-				newSolver.solveEquations(dir, new_x, new_b, nullptr, nullptr);
+				newSolver.solveEquations(dir, new_x, new_b);
 
 				for (size_t i = 0, size = new_x.size(); i < size; ++i) {
 					x.at(i) = new_x.at(i);
@@ -422,14 +422,8 @@ namespace storm {
 		}
         
         template<typename ValueType>
-        void TopologicalMinMaxLinearEquationSolver<ValueType>::multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
-            
-            // If scratch memory was not provided, we need to create it.
-            bool multiplyResultMemoryProvided = true;
-            if (multiplyResult == nullptr) {
-                multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
-                multiplyResultMemoryProvided = false;
-            }
+        void TopologicalMinMaxLinearEquationSolver<ValueType>::repeatedMultiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const {
+            std::unique_ptr<std::vector<ValueType>> multiplyResult = std::make_unique<std::vector<ValueType>>(this->A.getRowCount());
             
             // Now perform matrix-vector multiplication as long as we meet the bound of the formula.
             for (uint_fast64_t i = 0; i < n; ++i) {
@@ -445,10 +439,6 @@ namespace storm {
                 storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices());
                 
             }
-            
-            if (!multiplyResultMemoryProvided) {
-                delete multiplyResult;
-            }
         }
 
         template<typename ValueType>
diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.h b/src/solver/TopologicalMinMaxLinearEquationSolver.h
index e6cc87997..5eac4a976 100644
--- a/src/solver/TopologicalMinMaxLinearEquationSolver.h
+++ b/src/solver/TopologicalMinMaxLinearEquationSolver.h
@@ -32,9 +32,9 @@ namespace storm {
              */
             TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision = 1e-6, uint_fast64_t maximalNumberOfIterations = 20000, bool relative = true);
                         
-            virtual void solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
+            virtual void solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
             
-            virtual void multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const override;
+            virtual void repeatedMultiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const override;
 
         private:
             storm::storage::SparseMatrix<ValueType> const& A;
diff --git a/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp b/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp
index 4f86c25b1..94f49c2bf 100644
--- a/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp
+++ b/test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp
@@ -78,23 +78,23 @@ TEST(GmmxxMinMaxLinearEquationSolver, MatrixVectorMultiplication) {
     auto factory = storm::solver::GmmxxMinMaxLinearEquationSolverFactory<double>();
     auto solver = factory.create(A);
     
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 1));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 1));
     ASSERT_LT(std::abs(x[0] - 0.099), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 2));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 2));
     ASSERT_LT(std::abs(x[0] - 0.1881), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 20));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 20));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Maximize, x, nullptr, 1));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Maximize, x, nullptr, 1));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Maximize, x, nullptr, 20));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Maximize, x, nullptr, 20));
     ASSERT_LT(std::abs(x[0] - 0.9238082658), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
 }
 
diff --git a/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp b/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp
index fdf595193..4a61de149 100644
--- a/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp
+++ b/test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp
@@ -49,23 +49,23 @@ TEST(NativeMinMaxLinearEquationSolver, MatrixVectorMultiplication) {
     auto factory = storm::solver::NativeMinMaxLinearEquationSolverFactory<double>();
     auto solver = factory.create(A);
     
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 1));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 1));
     ASSERT_LT(std::abs(x[0] - 0.099), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 2));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 2));
     ASSERT_LT(std::abs(x[0] - 0.1881), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 20));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 20));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Maximize, x, nullptr, 1));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Maximize, x, nullptr, 1));
     ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
     
     x = {0, 1, 0};
-    ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Maximize, x, nullptr, 20));
+    ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Maximize, x, nullptr, 20));
     ASSERT_LT(std::abs(x[0] - 0.9238082658), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
 }