From be5fdeb636ebc049242e5d2c63446e8f4c653588 Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 28 Jun 2016 15:54:45 +0200 Subject: [PATCH] started working on internal auxiliary storage of solvers Former-commit-id: d895041c50a04955b365c0a192f8fd10adc44dfb --- .../csl/helper/SparseCtmcCslHelper.cpp | 9 +- .../helper/SparseMarkovAutomatonCslHelper.cpp | 7 +- .../prctl/helper/HybridMdpPrctlHelper.cpp | 6 +- .../prctl/helper/SparseDtmcPrctlHelper.cpp | 2 +- .../prctl/helper/SparseMdpPrctlHelper.cpp | 12 +- src/solver/EigenLinearEquationSolver.cpp | 31 +++- src/solver/EigenLinearEquationSolver.h | 12 +- .../EliminationLinearEquationSolver.cpp | 42 +++-- src/solver/EliminationLinearEquationSolver.h | 16 +- src/solver/GmmxxLinearEquationSolver.cpp | 121 +++++++++--- src/solver/GmmxxLinearEquationSolver.h | 26 ++- src/solver/LinearEquationSolver.cpp | 72 ++++++-- src/solver/LinearEquationSolver.h | 72 +++++++- src/solver/MinMaxLinearEquationSolver.cpp | 23 ++- src/solver/MinMaxLinearEquationSolver.h | 45 +++-- src/solver/NativeLinearEquationSolver.cpp | 153 ++++++++++------ src/solver/NativeLinearEquationSolver.h | 24 ++- .../StandardMinMaxLinearEquationSolver.cpp | 172 ++++++++++++------ .../StandardMinMaxLinearEquationSolver.h | 20 +- 19 files changed, 634 insertions(+), 231 deletions(-) diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index 4ccf61e64..186762a3d 100644 --- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -627,19 +627,19 @@ namespace storm { result = std::vector(values.size()); } } - std::vector multiplicationResult(result.size()); std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(uniformizedMatrix)); + solver->allocateAuxStorage(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 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 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 storm::storage::SparseMatrix SparseCtmcCslHelper::computeProbabilityMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector 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..b77a19825 100644 --- a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -86,6 +86,7 @@ namespace storm { } std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(aProbabilistic); + solver->allocateAuxStorage(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 markovianNonGoalValuesSwap(markovianNonGoalValues); - std::vector multiplicationResultScratchMemory(aProbabilistic.getRowCount()); - std::vector 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 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, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); std::unique_ptr> 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(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.template toAdd(), maybeStates, odd, x)); @@ -166,7 +166,7 @@ namespace storm { // Perform the matrix-vector multiplication. std::unique_ptr> 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(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), odd, x)); @@ -195,7 +195,7 @@ namespace storm { // Perform the matrix-vector multiplication. std::unique_ptr> 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(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), 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> 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..483169789 100644 --- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -49,7 +49,7 @@ namespace storm { std::vector subresult(maybeStates.getNumberOfSetBits()); std::unique_ptr> 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()); std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix); - solver->multiply(dir, result); + solver->repeatedMultiply(dir, result); return result; } @@ -211,7 +211,7 @@ namespace storm { std::vector result(rewardModel.getStateRewardVector()); std::unique_ptr> 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> 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() - 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() - 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 EigenLinearEquationSolver::EigenLinearEquationSolver(storm::storage::SparseMatrix&& A, EigenLinearEquationSolverSettings const& settings) : settings(settings) { + this->setMatrix(std::move(A)); + } + + template + void EigenLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix(A); + } + + template + void EigenLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + // Take ownership of the matrix so it is destroyed after we have translated it to Eigen's format. storm::storage::SparseMatrix localA(std::move(A)); - eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix(localA); + this->setMatrix(localA); } template - void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { // Map the input vectors to Eigen's format. auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); auto eigenB = Eigen::Matrix::Map(b.data(), b.size()); @@ -197,7 +208,7 @@ namespace storm { } template - void EigenLinearEquationSolver::multiply(std::vector& x, std::vector& result, std::vector const* b) const { + void EigenLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { // Typedef the map-type so we don't have to spell it out. typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; @@ -234,10 +245,20 @@ namespace storm { return settings; } + template + uint64_t EigenLinearEquationSolver::getMatrixRowCount() const { + return eigenA->rows(); + } + + template + uint64_t EigenLinearEquationSolver::getMatrixColumnCount() const { + return eigenA->cols(); + } + // Specialization form storm::RationalNumber template<> - void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { // Map the input vectors to Eigen's format. auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); auto eigenB = Eigen::Matrix::Map(b.data(), b.size()); @@ -250,7 +271,7 @@ namespace storm { // Specialization form storm::RationalFunction template<> - void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { // Map the input vectors to Eigen's format. auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); auto eigenB = Eigen::Matrix::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 const& A, EigenLinearEquationSolverSettings const& settings = EigenLinearEquationSolverSettings()); EigenLinearEquationSolver(storm::storage::SparseMatrix&& A, EigenLinearEquationSolverSettings const& settings = EigenLinearEquationSolverSettings()); - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual void solveEquations(std::vector& x, std::vector const& b) const override; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; EigenLinearEquationSolverSettings& getSettings(); EigenLinearEquationSolverSettings 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> 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 - EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A, EliminationLinearEquationSolverSettings const& settings) : localA(nullptr), A(A), settings(settings) { - // Intentionally left empty. + EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A, EliminationLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings) { + this->setMatrix(A); } template - EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix&& A, EliminationLinearEquationSolverSettings const& settings) : localA(std::make_unique>(std::move(A))), A(*localA), settings(settings) { - // Intentionally left empty. + EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix&& A, EliminationLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings) { + this->setMatrix(std::move(A)); } template - void EliminationLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { - STORM_LOG_WARN_COND(multiplyResult == nullptr, "Providing scratch memory will not improve the performance of this solver."); - + void EliminationLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + this->A = &A; + localA.reset(); + } + + template + void EliminationLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + localA = std::make_unique>(std::move(A)); + this->A = localA.get(); + } + + template + void EliminationLinearEquationSolver::solveEquations(std::vector& x, std::vector 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 - void EliminationLinearEquationSolver::multiply(std::vector& x, std::vector& result, std::vector const* b) const { + void EliminationLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& 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 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 + uint64_t EliminationLinearEquationSolver::getMatrixRowCount() const { + return this->A->getRowCount(); + } + + template + uint64_t EliminationLinearEquationSolver::getMatrixColumnCount() const { + return this->A->getColumnCount(); + } + template std::unique_ptr> EliminationLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::make_unique>(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 const& A, EliminationLinearEquationSolverSettings const& settings = EliminationLinearEquationSolverSettings()); EliminationLinearEquationSolver(storm::storage::SparseMatrix&& A, EliminationLinearEquationSolverSettings const& settings = EliminationLinearEquationSolverSettings()); - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual void solveEquations(std::vector& x, std::vector const& b) const override; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; EliminationLinearEquationSolverSettings& getSettings(); EliminationLinearEquationSolverSettings 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> 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 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 const* A; // The settings used by the solver. EliminationLinearEquationSolverSettings settings; diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index d29cda389..139a969dd 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -111,17 +111,31 @@ namespace storm { } template - GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(A), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A)), settings(settings) { - // Intentionally left empty. + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiStorage(nullptr) { + this->setMatrix(A); } template - GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings) : localA(std::make_unique>(std::move(A))), A(*localA), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(*localA)), settings(settings) { - // Intentionally left empty. + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiStorage(nullptr) { + this->setMatrix(std::move(A)); } template - void GmmxxLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void GmmxxLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + localA.reset(); + this->A = &A; + gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + } + + template + void GmmxxLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + localA = std::make_unique>(std::move(A)); + this->A = localA.get(); + gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(*localA); + } + + template + void GmmxxLinearEquationSolver::solveEquations(std::vector& x, std::vector 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::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 - void GmmxxLinearEquationSolver::multiply(std::vector& x, std::vector& result, std::vector const* b) const { + void GmmxxLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { if (b) { gmm::mult_add(*gmmxxMatrix, x, *b, result); } else { @@ -187,25 +201,20 @@ namespace storm { } template - uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { + bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliaryJacobiStorage = std::make_unique>(x.size()); + } + // Get a Jacobi decomposition of the matrix A. std::pair, std::vector> jacobiDecomposition = A.getJacobiDecomposition(); // Convert the LU matrix to gmm++'s format. std::unique_ptr> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(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* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; std::vector* currentX = &x; - - // Target vector for precision calculation. - std::vector tmpX(x.size()); + std::vector* nextX = auxiliaryJacobiStorage.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()), tmpX); - storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX); + gmm::mult_add(*gmmLU, gmm::scaled(*currentX, -storm::utility::one()), 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 == auxiliaryJacobiStorage.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 storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliaryJacobiStorage.reset(); } return iterationCount; @@ -251,6 +259,67 @@ namespace storm { return settings; } + template + bool GmmxxLinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi) { + if (!auxiliaryJacobiStorage) { + auxiliaryJacobiStorage = std::make_unique>(this->getMatrixRowCount()); + result = true; + } + } + } + result |= LinearEquationSolver::allocateAuxStorage(operation); + return result; + } + + template + bool GmmxxLinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliaryJacobiStorage) { + result = true; + } + auxiliaryJacobiStorage.reset(); + } + result |= LinearEquationSolver::deallocateAuxStorage(operation); + return result; + } + + template + bool GmmxxLinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliaryJacobiStorage != nullptr) { + result = auxiliaryJacobiStorage->size() != this->getMatrixColumnCount(); + auxiliaryJacobiStorage->resize(this->getMatrixRowCount()); + } + } + result |= LinearEquationSolver::reallocateAuxStorage(operation); + return result; + } + + template + bool GmmxxLinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + result |= auxiliaryJacobiStorage != nullptr; + } + result |= LinearEquationSolver::hasAuxStorage(operation); + return result; + } + + template + uint64_t GmmxxLinearEquationSolver::getMatrixRowCount() const { + return this->A->getRowCount(); + } + + template + uint64_t GmmxxLinearEquationSolver::getMatrixColumnCount() const { + return this->A->getColumnCount(); + } + template std::unique_ptr> GmmxxLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::make_unique>(matrix, settings); diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 4425ec940..57386677d 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -89,12 +89,20 @@ namespace storm { GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings = GmmxxLinearEquationSolverSettings()); GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings = GmmxxLinearEquationSolverSettings()); - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual void solveEquations(std::vector& x, std::vector const& b) const override; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; GmmxxLinearEquationSolverSettings& getSettings(); GmmxxLinearEquationSolverSettings const& getSettings() const; + virtual bool allocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool hasAuxStorage(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 const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const; + uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector 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> 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 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 const* A; // The (gmm++) matrix associated with this equation solver. std::unique_ptr> gmmxxMatrix; // The settings used by the solver. GmmxxLinearEquationSolverSettings settings; + + // Auxiliary storage for the Jacobi method. + mutable std::unique_ptr> auxiliaryJacobiStorage; }; template diff --git a/src/solver/LinearEquationSolver.cpp b/src/solver/LinearEquationSolver.cpp index 84f5f84e5..feafdfde2 100644 --- a/src/solver/LinearEquationSolver.cpp +++ b/src/solver/LinearEquationSolver.cpp @@ -14,36 +14,80 @@ namespace storm { namespace solver { template - void LinearEquationSolver::repeatedMultiply(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { + LinearEquationSolver::LinearEquationSolver() : auxiliaryRepeatedMultiplyStorage(nullptr) { + // Intentionally left empty. + } + + template + void LinearEquationSolver::repeatedMultiply(std::vector& x, std::vector const* b, uint_fast64_t n) const { + bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::MultiplyRepeatedly); + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage = std::make_unique>(x.size()); + } // Set up some temporary variables so that we can just swap pointers instead of copying the result after // each iteration. std::vector* currentX = &x; - - bool multiplyResultProvided = true; - std::vector* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; + std::vector* nextX = auxiliaryRepeatedMultiplyStorage.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 == auxiliaryRepeatedMultiplyStorage.get()) { std::swap(x, *currentX); } + + // If we allocated auxiliary storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage.reset(); + } + } + + template + bool LinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { + if (this->hasAuxStorage(operation)) { + return false; + } - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; + auxiliaryRepeatedMultiplyStorage = std::make_unique>(this->getMatrixColumnCount()); + return true; + } + + template + bool LinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { + bool result = auxiliaryRepeatedMultiplyStorage != nullptr; + if (result) { + auxiliaryRepeatedMultiplyStorage.reset(); + } + return result; + } + return false; + } + + template + bool LinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { + if (auxiliaryRepeatedMultiplyStorage != nullptr) { + result = auxiliaryRepeatedMultiplyStorage->size() != this->getMatrixColumnCount(); + auxiliaryRepeatedMultiplyStorage->resize(this->getMatrixColumnCount()); + } + } + return result; + } + + template + bool LinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { + return auxiliaryRepeatedMultiplyStorage != nullptr; } + return false; } template diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h index 9632ddb2b..2e18ecd0c 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 LinearEquationSolver : public AbstractEquationSolver { public: + LinearEquationSolver(); + virtual ~LinearEquationSolver() { // Intentionally left empty. } + virtual void setMatrix(storm::storage::SparseMatrix const& A) = 0; + virtual void setMatrix(storm::storage::SparseMatrix&& 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& x, std::vector const& b, std::vector* multiplyResult = nullptr) const = 0; + virtual void solveEquations(std::vector& x, std::vector 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& x, std::vector& result, std::vector const* b = nullptr) const = 0; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& 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& x, std::vector const* b, uint_fast64_t n) const; + + // 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 allocateAuxStorage(LinearEquationSolverOperation operation); + + /*! + * Destroys previously allocated auxiliary storage for the provided operation. + * + * @return True iff auxiliary storage was deallocated. */ - void repeatedMultiply(std::vector& x, std::vector const* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const; + virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation); + + /*! + * If the matrix dimensions changed and auxiliary storage was allocated, this function needs to be called to + * update the auxiliary storage. + * + * @return True iff the auxiliary storage was reallocated. + */ + virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation); + + /*! + * 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 hasAuxStorage(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 storage for repeated matrix-vector multiplication. + mutable std::unique_ptr> auxiliaryRepeatedMultiplyStorage; }; template diff --git a/src/solver/MinMaxLinearEquationSolver.cpp b/src/solver/MinMaxLinearEquationSolver.cpp index bd522834f..bb070c036 100644 --- a/src/solver/MinMaxLinearEquationSolver.cpp +++ b/src/solver/MinMaxLinearEquationSolver.cpp @@ -28,15 +28,15 @@ namespace storm { } template - void MinMaxLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void MinMaxLinearEquationSolver::solveEquations(std::vector& x, std::vector 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 - void MinMaxLinearEquationSolver::multiply( std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const { + void MinMaxLinearEquationSolver::repeatedMultiply( std::vector& x, std::vector* 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 @@ -79,6 +79,21 @@ namespace storm { return std::move(scheduler.get()); } + template + bool MinMaxLinearEquationSolver::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + return false; + } + + template + bool MinMaxLinearEquationSolver::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + return false; + } + + template + bool MinMaxLinearEquationSolver::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + return false; + } + template MinMaxLinearEquationSolverFactory::MinMaxLinearEquationSolverFactory(bool trackScheduler) : trackScheduler(trackScheduler) { // Intentionally left empty. diff --git a/src/solver/MinMaxLinearEquationSolver.h b/src/solver/MinMaxLinearEquationSolver.h index 358ef0be1..565306756 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& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const = 0; + virtual void solveEquations(OptimizationDirection d, std::vector& x, std::vector const& b) const = 0; /*! * Behaves the same as the other variant of solveEquations, 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& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const; + void solveEquations(std::vector& x, std::vector 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& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const = 0; + virtual void repeatedMultiply(OptimizationDirection d, std::vector& x, std::vector* b, uint_fast64_t n = 1) const = 0; /*! * Behaves the same as the other variant of multiply, 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& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const; + virtual void repeatedMultiply(std::vector& x, std::vector* 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 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 allocateAuxStorage(MinMaxLinearEquationSolverOperation operation); + + /*! + * Destroys previously allocated auxiliary storage for the provided operation. + * + * @return True iff auxiliary storage was deallocated. + */ + virtual bool deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation); + + /*! + * 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 hasAuxStorage(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..0eb5d758e 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -84,70 +84,62 @@ namespace storm { } template - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(A), settings(settings) { - // Intentionally left empty. + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingStorage(nullptr) { + this->setMatrix(A); } template - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings) : localA(std::make_unique>(std::move(A))), A(*localA), settings(settings) { - // Intentionally left empty. + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingStorage(nullptr) { + this->setMatrix(std::move(A)); } template - void NativeLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void NativeLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + localA.reset(); + this->A = &A; + } + + template + void NativeLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + localA = std::make_unique>(std::move(A)); + this->A = localA.get(); + } + + template + void NativeLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { + bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliarySolvingStorage = std::make_unique>(x.size()); + } + if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::GaussSeidel) { // Define the omega used for SOR. ValueType omega = this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR ? this->getSettings().getOmega() : storm::utility::one(); - // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with. - bool tmpXProvided = true; - std::vector* tmpX = multiplyResult; - if (multiplyResult == nullptr) { - tmpX = new std::vector(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(x, *tmpX, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)); + converged = storm::utility::vector::equalModuloPrecision(*auxiliarySolvingStorage, x, static_cast(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; + *auxiliarySolvingStorage = 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, std::vector> jacobiDecomposition = A.getJacobiDecomposition(); + std::pair, std::vector> 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* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; std::vector* currentX = &x; - - // Target vector for multiplication. - std::vector tmpX(x.size()); + std::vector* nextX = auxiliarySolvingStorage.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(*currentX, *nextX, static_cast(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 == auxiliarySolvingStorage.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 storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliarySolvingStorage.reset(); } } template - void NativeLinearEquationSolver::multiply(std::vector& x, std::vector& result, std::vector const* b) const { + void NativeLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& 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 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 + bool NativeLinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (!auxiliarySolvingStorage) { + auxiliarySolvingStorage = std::make_unique>(this->getMatrixRowCount()); + result = true; + } + } + result |= LinearEquationSolver::allocateAuxStorage(operation); + return result; + } + + template + bool NativeLinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliarySolvingStorage) { + result = true; + } + auxiliarySolvingStorage.reset(); + } + result |= LinearEquationSolver::deallocateAuxStorage(operation); + return result; + } + + template + bool NativeLinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliarySolvingStorage) { + result = auxiliarySolvingStorage->size() != this->getMatrixColumnCount(); + auxiliarySolvingStorage->resize(this->getMatrixRowCount()); + } + } + result |= LinearEquationSolver::reallocateAuxStorage(operation); + return result; + } + + template + bool NativeLinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + result |= auxiliarySolvingStorage != nullptr; + } + result |= LinearEquationSolver::hasAuxStorage(operation); + return result; + } + + template + uint64_t NativeLinearEquationSolver::getMatrixRowCount() const { + return this->A->getRowCount(); + } + + template + uint64_t NativeLinearEquationSolver::getMatrixColumnCount() const { + return this->A->getColumnCount(); + } + template std::unique_ptr> NativeLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::make_unique>(matrix, settings); diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index 235985df6..19722debe 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -46,23 +46,37 @@ namespace storm { NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings = NativeLinearEquationSolverSettings()); NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings = NativeLinearEquationSolverSettings()); - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual void solveEquations(std::vector& x, std::vector const& b) const override; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; NativeLinearEquationSolverSettings& getSettings(); NativeLinearEquationSolverSettings const& getSettings() const; + virtual bool allocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool hasAuxStorage(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> 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 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 const* A; // The settings used by the solver. NativeLinearEquationSolverSettings settings; + + // Auxiliary storage for the equation solving methods. + mutable std::unique_ptr> auxiliarySolvingStorage; }; template diff --git a/src/solver/StandardMinMaxLinearEquationSolver.cpp b/src/solver/StandardMinMaxLinearEquationSolver.cpp index db507ef68..2af6527e7 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/solver/StandardMinMaxLinearEquationSolver.cpp @@ -84,45 +84,38 @@ namespace storm { } template - void StandardMinMaxLinearEquationSolver::solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void StandardMinMaxLinearEquationSolver::solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const { switch (this->getSettings().getSolutionMethod()) { - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: solveEquationsValueIteration(dir, x, b, multiplyResult, newX); break; - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: solveEquationsPolicyIteration(dir, x, b, multiplyResult, newX); break; + case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: solveEquationsValueIteration(dir, x, b); break; + case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: solveEquationsPolicyIteration(dir, x, b); break; } } template - void StandardMinMaxLinearEquationSolver::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { - - // Create scratch memory if none was provided. - bool multiplyResultMemoryProvided = multiplyResult != nullptr; - if (multiplyResult == nullptr) { - multiplyResult = new std::vector(this->A.getRowCount()); - } - + void StandardMinMaxLinearEquationSolver::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { // Create the initial scheduler. std::vector scheduler(this->A.getRowGroupCount()); // Create a vector for storing the right-hand side of the inner equation system. std::vector subB(this->A.getRowGroupCount()); - - // Create a vector that the inner equation solver can use as scratch memory. - std::vector deterministicMultiplyResult(this->A.getRowGroupCount()); + + // Resolve the nondeterminism according to the current scheduler. + storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(scheduler, true); + submatrix.convertToEquationSystem(); + storm::utility::vector::selectVectorValues(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->allocateAuxStorage(LinearEquationSolverOperation::SolveEquations); Status status = Status::InProgress; uint64_t iterations = 0; do { - // Resolve the nondeterminism according to the current scheduler. - storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(scheduler, true); - submatrix.convertToEquationSystem(); - storm::utility::vector::selectVectorValues(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(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(std::move(scheduler)); } - - if (!multiplyResultMemoryProvided) { - delete multiplyResult; - } } template @@ -186,22 +181,16 @@ namespace storm { } template - void StandardMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void StandardMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { std::unique_ptr> solver = linearEquationSolverFactory->create(A); - - // Create scratch memory if none was provided. - bool multiplyResultMemoryProvided = multiplyResult != nullptr; - if (multiplyResult == nullptr) { - multiplyResult = new std::vector(this->A.getRowCount()); - } - std::vector* currentX = &x; - bool xMemoryProvided = newX != nullptr; - if (newX == nullptr) { - newX = new std::vector(x.size()); + bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliarySolvingMultiplyStorage = std::make_unique>(this->A.getRowCount()); + auxiliarySolvingVectorStorage = std::make_unique>(x.size()); } - // Keep track of which of the vectors for x is the auxiliary copy. - std::vector* copyX = newX; + std::vector* currentX = &x; + std::vector* newX = auxiliarySolvingVectorStorage.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 +198,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, *auxiliarySolvingMultiplyStorage); // 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, *auxiliarySolvingMultiplyStorage, *newX, this->A.getRowGroupIndices()); // Determine whether the method converged. if (storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) { @@ -229,39 +218,37 @@ 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 == auxiliarySolvingVectorStorage.get()) { std::swap(x, *currentX); } - // Dispose of allocated scratch memory. - if (!xMemoryProvided) { - delete copyX; - } - if (!multiplyResultMemoryProvided) { - delete multiplyResult; + // If we allocated auxiliary storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliarySolvingMultiplyStorage.reset(); + auxiliarySolvingVectorStorage.reset(); } } template - void StandardMinMaxLinearEquationSolver::multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const { - std::unique_ptr> solver = linearEquationSolverFactory->create(A); - - // If scratch memory was not provided, we create it. - bool multiplyResultMemoryProvided = multiplyResult != nullptr; - if (!multiplyResult) { - multiplyResult = new std::vector(this->A.getRowCount()); + void StandardMinMaxLinearEquationSolver::multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const { + bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly); + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage = std::make_unique>(this->A.getRowCount()); } + std::unique_ptr> solver = linearEquationSolverFactory->create(A); + for (uint64_t i = 0; i < n; ++i) { - solver->multiply(x, *multiplyResult, b); + solver->multiply(x, b, *auxiliaryRepeatedMultiplyStorage); // 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, *auxiliaryRepeatedMultiplyStorage, x, this->A.getRowGroupIndices()); } - if (!multiplyResultMemoryProvided) { - delete multiplyResult; + // If we allocated auxiliary storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage.reset(); } } @@ -298,6 +285,75 @@ namespace storm { return settings; } + template + bool StandardMinMaxLinearEquationSolver::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool result = false; + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { + if (!auxiliarySolvingMultiplyStorage) { + result = true; + auxiliarySolvingMultiplyStorage = std::make_unique>(this->A.getRowCount()); + } + if (!auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage = std::make_unique>(this->A.getRowGroupCount()); + } + } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration) { + if (!auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage = std::make_unique>(this->A.getRowGroupCount()); + } + } else { + STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method."); + } + } else { + if (!auxiliaryRepeatedMultiplyStorage) { + result = true; + auxiliaryRepeatedMultiplyStorage = std::make_unique>(this->A.getRowCount()); + } + } + return result; + } + + template + bool StandardMinMaxLinearEquationSolver::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool result = false; + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { + if (auxiliarySolvingMultiplyStorage) { + result = true; + auxiliarySolvingMultiplyStorage.reset(); + } + if (auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage.reset(); + } + } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration) { + // Nothing to do in this case. + } else { + STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method."); + } + } else { + if (auxiliaryRepeatedMultiplyStorage) { + result = true; + auxiliaryRepeatedMultiplyStorage.reset(); + } + } + } + + template + bool StandardMinMaxLinearEquationSolver::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { + return auxiliarySolvingMultiplyStorage && auxiliarySolvingVectorStorage; + } else { + return false; + } + } else { + return static_cast(auxiliaryRepeatedMultiplyStorage); + } + } + template StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler), linearEquationSolverFactory(nullptr) { // Intentionally left empty. diff --git a/src/solver/StandardMinMaxLinearEquationSolver.h b/src/solver/StandardMinMaxLinearEquationSolver.h index 36b2aa458..523266c04 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/solver/StandardMinMaxLinearEquationSolver.h @@ -38,15 +38,19 @@ namespace storm { StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); - virtual void solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const override; - virtual void multiply(OptimizationDirection dir, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; + virtual void solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const override; + virtual void multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const override; StandardMinMaxLinearEquationSolverSettings const& getSettings() const; StandardMinMaxLinearEquationSolverSettings& getSettings(); + virtual bool allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; + virtual bool deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; + virtual bool hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const override; + private: - void solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const; - void solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const; + void solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + void solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const; @@ -70,6 +74,14 @@ 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 const& A; + + // Auxiliary storage for equation solving. + // Auxiliary storage for repeated matrix-vector multiplication. + mutable std::unique_ptr> auxiliarySolvingMultiplyStorage; + mutable std::unique_ptr> auxiliarySolvingVectorStorage; + + // Auxiliary storage for repeated matrix-vector multiplication. + mutable std::unique_ptr> auxiliaryRepeatedMultiplyStorage; }; template