From be5fdeb636ebc049242e5d2c63446e8f4c653588 Mon Sep 17 00:00:00 2001 From: dehnert <dehnert@cs.rwth-aachen.de> Date: Tue, 28 Jun 2016 15:54:45 +0200 Subject: [PATCH 1/2] 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<ValueType>(values.size()); } } - std::vector<ValueType> multiplicationResult(result.size()); std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> 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<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..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<storm::solver::MinMaxLinearEquationSolver<ValueType>> 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<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..483169789 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); 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..139a969dd 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), auxiliaryJacobiStorage(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), auxiliaryJacobiStorage(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 allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliaryJacobiStorage = std::make_unique<std::vector<ValueType>>(x.size()); + } + // 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 = 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<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 == 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<typename ValueType> + bool GmmxxLinearEquationSolver<ValueType>::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) { + if (!auxiliaryJacobiStorage) { + auxiliaryJacobiStorage = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount()); + result = true; + } + } + } + result |= LinearEquationSolver<ValueType>::allocateAuxStorage(operation); + return result; + } + + template<typename ValueType> + bool GmmxxLinearEquationSolver<ValueType>::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliaryJacobiStorage) { + result = true; + } + auxiliaryJacobiStorage.reset(); + } + result |= LinearEquationSolver<ValueType>::deallocateAuxStorage(operation); + return result; + } + + template<typename ValueType> + bool GmmxxLinearEquationSolver<ValueType>::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliaryJacobiStorage != nullptr) { + result = auxiliaryJacobiStorage->size() != this->getMatrixColumnCount(); + auxiliaryJacobiStorage->resize(this->getMatrixRowCount()); + } + } + result |= LinearEquationSolver<ValueType>::reallocateAuxStorage(operation); + return result; + } + + template<typename ValueType> + bool GmmxxLinearEquationSolver<ValueType>::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + result |= auxiliaryJacobiStorage != nullptr; + } + result |= LinearEquationSolver<ValueType>::hasAuxStorage(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..57386677d 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 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<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>> auxiliaryJacobiStorage; }; template<typename ValueType> 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<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() : auxiliaryRepeatedMultiplyStorage(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 allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::MultiplyRepeatedly); + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage = std::make_unique<std::vector<ValueType>>(x.size()); + } // 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 = 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<typename ValueType> + bool LinearEquationSolver<ValueType>::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<std::vector<ValueType>>(this->getMatrixColumnCount()); + return true; + } + + template<typename ValueType> + bool LinearEquationSolver<ValueType>::deallocateAuxStorage(LinearEquationSolverOperation operation) { + if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { + bool result = auxiliaryRepeatedMultiplyStorage != nullptr; + if (result) { + auxiliaryRepeatedMultiplyStorage.reset(); + } + return result; + } + return false; + } + + template<typename ValueType> + bool LinearEquationSolver<ValueType>::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<typename ValueType> + bool LinearEquationSolver<ValueType>::hasAuxStorage(LinearEquationSolverOperation operation) const { + if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { + return auxiliaryRepeatedMultiplyStorage != nullptr; } + return false; } template<typename ValueType> 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 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 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<ValueType>& x, std::vector<ValueType> const* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* 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<std::vector<ValueType>> auxiliaryRepeatedMultiplyStorage; }; template<typename ValueType> 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<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>::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + return false; + } + + template<typename ValueType> + bool MinMaxLinearEquationSolver<ValueType>::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + return false; + } + + template<typename ValueType> + bool MinMaxLinearEquationSolver<ValueType>::hasAuxStorage(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..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<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 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<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), auxiliarySolvingStorage(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), auxiliarySolvingStorage(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->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliarySolvingStorage = std::make_unique<std::vector<ValueType>>(x.size()); + } + 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>(*auxiliarySolvingStorage, 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; + *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<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 = 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<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 == 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<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>::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (!auxiliarySolvingStorage) { + auxiliarySolvingStorage = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount()); + result = true; + } + } + result |= LinearEquationSolver<ValueType>::allocateAuxStorage(operation); + return result; + } + + template<typename ValueType> + bool NativeLinearEquationSolver<ValueType>::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliarySolvingStorage) { + result = true; + } + auxiliarySolvingStorage.reset(); + } + result |= LinearEquationSolver<ValueType>::deallocateAuxStorage(operation); + return result; + } + + template<typename ValueType> + bool NativeLinearEquationSolver<ValueType>::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliarySolvingStorage) { + result = auxiliarySolvingStorage->size() != this->getMatrixColumnCount(); + auxiliarySolvingStorage->resize(this->getMatrixRowCount()); + } + } + result |= LinearEquationSolver<ValueType>::reallocateAuxStorage(operation); + return result; + } + + template<typename ValueType> + bool NativeLinearEquationSolver<ValueType>::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + result |= auxiliarySolvingStorage != nullptr; + } + result |= LinearEquationSolver<ValueType>::hasAuxStorage(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..19722debe 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 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<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 storage for the equation solving methods. + mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingStorage; }; template<typename ValueType> 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<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->allocateAuxStorage(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,16 @@ 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 allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliarySolvingMultiplyStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); + auxiliarySolvingVectorStorage = std::make_unique<std::vector<ValueType>>(x.size()); } - // 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 = 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<ValueType>(*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<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>::multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const { + bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly); + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); } + 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, *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<typename ValueType> + bool StandardMinMaxLinearEquationSolver<ValueType>::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool result = false; + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) { + if (!auxiliarySolvingMultiplyStorage) { + result = true; + auxiliarySolvingMultiplyStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); + } + if (!auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowGroupCount()); + } + } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) { + if (!auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowGroupCount()); + } + } else { + STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method."); + } + } else { + if (!auxiliaryRepeatedMultiplyStorage) { + result = true; + auxiliaryRepeatedMultiplyStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); + } + } + return result; + } + + template<typename ValueType> + bool StandardMinMaxLinearEquationSolver<ValueType>::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool result = false; + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) { + if (auxiliarySolvingMultiplyStorage) { + result = true; + auxiliarySolvingMultiplyStorage.reset(); + } + if (auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage.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 (auxiliaryRepeatedMultiplyStorage) { + result = true; + auxiliaryRepeatedMultiplyStorage.reset(); + } + } + } + + template<typename ValueType> + bool StandardMinMaxLinearEquationSolver<ValueType>::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) { + return auxiliarySolvingMultiplyStorage && auxiliarySolvingVectorStorage; + } else { + return false; + } + } else { + return static_cast<bool>(auxiliaryRepeatedMultiplyStorage); + } + } + 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..523266c04 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 multiply(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 allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; + virtual bool deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; + virtual bool hasAuxStorage(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,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<ValueType> const& A; + + // Auxiliary storage for equation solving. + // Auxiliary storage for repeated matrix-vector multiplication. + mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingMultiplyStorage; + mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingVectorStorage; + + // Auxiliary storage for repeated matrix-vector multiplication. + mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRepeatedMultiplyStorage; }; template<typename ValueType> From 83c4b1647c7a9918926a82440e2eda551484813b Mon Sep 17 00:00:00 2001 From: dehnert <dehnert@cs.rwth-aachen.de> Date: Tue, 28 Jun 2016 18:48:14 +0200 Subject: [PATCH 2/2] solvers now can allocated auxiliary memory Former-commit-id: 76dc1a16798b92567a2ce6c4f9cae54ca9343677 --- .../csl/helper/SparseCtmcCslHelper.cpp | 2 +- .../helper/SparseMarkovAutomatonCslHelper.cpp | 2 +- .../prctl/helper/SparseMdpPrctlHelper.cpp | 2 +- src/solver/GmmxxLinearEquationSolver.cpp | 52 ++++++------ src/solver/GmmxxLinearEquationSolver.h | 10 +-- src/solver/LinearEquationSolver.cpp | 50 ++++++----- src/solver/LinearEquationSolver.h | 36 ++++---- src/solver/MinMaxLinearEquationSolver.cpp | 6 +- src/solver/MinMaxLinearEquationSolver.h | 6 +- src/solver/NativeLinearEquationSolver.cpp | 52 ++++++------ src/solver/NativeLinearEquationSolver.h | 12 +-- .../StandardMinMaxLinearEquationSolver.cpp | 82 +++++++++---------- .../StandardMinMaxLinearEquationSolver.h | 19 ++--- .../TopologicalMinMaxLinearEquationSolver.cpp | 18 +--- .../TopologicalMinMaxLinearEquationSolver.h | 4 +- .../GmmxxMinMaxLinearEquationSolverTest.cpp | 10 +-- .../NativeMinMaxLinearEquationSolverTest.cpp | 10 +-- 17 files changed, 178 insertions(+), 195 deletions(-) diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index 186762a3d..dab1449de 100644 --- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -629,7 +629,7 @@ namespace storm { } std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(uniformizedMatrix)); - solver->allocateAuxStorage(storm::solver::LinearEquationSolverOperation::MultiplyRepeatedly); + solver->allocateAuxMemory(storm::solver::LinearEquationSolverOperation::MultiplyRepeatedly); if (!useMixedPoissonProbabilities && std::get<0>(foxGlynnResult) > 1) { // Perform the matrix-vector multiplications (without adding). diff --git a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index b77a19825..7817dc644 100644 --- a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -86,7 +86,7 @@ namespace storm { } std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(aProbabilistic); - solver->allocateAuxStorage(storm::solver::MinMaxLinearEquationSolverOperation::SolveEquations); + solver->allocateAuxMemory(storm::solver::MinMaxLinearEquationSolverOperation::SolveEquations); // Perform the actual value iteration // * loop until the step bound has been reached diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index 483169789..75fbfdfea 100644 --- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -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->repeatedMultiply(dir, result); + solver->repeatedMultiply(dir, result, nullptr, 1); return result; } diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index 139a969dd..76d38f726 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -111,12 +111,12 @@ namespace storm { } template<typename ValueType> - GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiStorage(nullptr) { + 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(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiStorage(nullptr) { + 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)); } @@ -202,9 +202,9 @@ 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) const { - bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); - if (allocatedAuxStorage) { - auxiliaryJacobiStorage = std::make_unique<std::vector<ValueType>>(x.size()); + bool allocatedAuxMemory = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations); + if (allocatedAuxMemory) { + this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations); } // Get a Jacobi decomposition of the matrix A. @@ -214,7 +214,7 @@ namespace storm { std::unique_ptr<gmm::csr_matrix<ValueType>> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.first)); std::vector<ValueType>* currentX = &x; - std::vector<ValueType>* nextX = auxiliaryJacobiStorage.get(); + std::vector<ValueType>* nextX = auxiliaryJacobiMemory.get(); // Set up additional environment variables. uint_fast64_t iterationCount = 0; @@ -237,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 == auxiliaryJacobiStorage.get()) { + if (currentX == auxiliaryJacobiMemory.get()) { std::swap(x, *currentX); } - // If we allocated auxiliary storage, we need to dispose of it now. - if (allocatedAuxStorage) { - auxiliaryJacobiStorage.reset(); + // If we allocated auxiliary memory, we need to dispose of it now. + if (allocatedAuxMemory) { + this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations); } return iterationCount; @@ -260,53 +260,53 @@ namespace storm { } template<typename ValueType> - bool GmmxxLinearEquationSolver<ValueType>::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool GmmxxLinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { if (this->getSettings().getSolutionMethod() == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) { - if (!auxiliaryJacobiStorage) { - auxiliaryJacobiStorage = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount()); + if (!auxiliaryJacobiMemory) { + auxiliaryJacobiMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount()); result = true; } } } - result |= LinearEquationSolver<ValueType>::allocateAuxStorage(operation); + result |= LinearEquationSolver<ValueType>::allocateAuxMemory(operation); return result; } template<typename ValueType> - bool GmmxxLinearEquationSolver<ValueType>::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool GmmxxLinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliaryJacobiStorage) { + if (auxiliaryJacobiMemory) { result = true; + auxiliaryJacobiMemory.reset(); } - auxiliaryJacobiStorage.reset(); } - result |= LinearEquationSolver<ValueType>::deallocateAuxStorage(operation); + result |= LinearEquationSolver<ValueType>::deallocateAuxMemory(operation); return result; } template<typename ValueType> - bool GmmxxLinearEquationSolver<ValueType>::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool GmmxxLinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliaryJacobiStorage != nullptr) { - result = auxiliaryJacobiStorage->size() != this->getMatrixColumnCount(); - auxiliaryJacobiStorage->resize(this->getMatrixRowCount()); + if (auxiliaryJacobiMemory) { + result = auxiliaryJacobiMemory->size() != this->getMatrixColumnCount(); + auxiliaryJacobiMemory->resize(this->getMatrixRowCount()); } } - result |= LinearEquationSolver<ValueType>::reallocateAuxStorage(operation); + result |= LinearEquationSolver<ValueType>::reallocateAuxMemory(operation); return result; } template<typename ValueType> - bool GmmxxLinearEquationSolver<ValueType>::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool GmmxxLinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - result |= auxiliaryJacobiStorage != nullptr; + result |= static_cast<bool>(auxiliaryJacobiMemory); } - result |= LinearEquationSolver<ValueType>::hasAuxStorage(operation); + result |= LinearEquationSolver<ValueType>::hasAuxMemory(operation); return result; } diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 57386677d..e7ab4c828 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -98,10 +98,10 @@ namespace storm { GmmxxLinearEquationSolverSettings<ValueType>& getSettings(); GmmxxLinearEquationSolverSettings<ValueType> 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; + 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: /*! @@ -136,7 +136,7 @@ namespace storm { GmmxxLinearEquationSolverSettings<ValueType> settings; // Auxiliary storage for the Jacobi method. - mutable std::unique_ptr<std::vector<ValueType>> auxiliaryJacobiStorage; + mutable std::unique_ptr<std::vector<ValueType>> auxiliaryJacobiMemory; }; template<typename ValueType> diff --git a/src/solver/LinearEquationSolver.cpp b/src/solver/LinearEquationSolver.cpp index feafdfde2..81dc4124b 100644 --- a/src/solver/LinearEquationSolver.cpp +++ b/src/solver/LinearEquationSolver.cpp @@ -14,21 +14,21 @@ namespace storm { namespace solver { template<typename ValueType> - LinearEquationSolver<ValueType>::LinearEquationSolver() : auxiliaryRepeatedMultiplyStorage(nullptr) { + 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 allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::MultiplyRepeatedly); - if (allocatedAuxStorage) { - auxiliaryRepeatedMultiplyStorage = std::make_unique<std::vector<ValueType>>(x.size()); + 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; - std::vector<ValueType>* nextX = auxiliaryRepeatedMultiplyStorage.get(); + 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) { @@ -38,54 +38,52 @@ namespace storm { // 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 == auxiliaryRepeatedMultiplyStorage.get()) { + if (currentX == auxiliaryRepeatedMultiplyMemory.get()) { std::swap(x, *currentX); } - // If we allocated auxiliary storage, we need to dispose of it now. - if (allocatedAuxStorage) { - auxiliaryRepeatedMultiplyStorage.reset(); + // If we allocated auxiliary memory, we need to dispose of it now. + if (allocatedAuxMemory) { + this->deallocateAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly); } } template<typename ValueType> - bool LinearEquationSolver<ValueType>::allocateAuxStorage(LinearEquationSolverOperation operation) { - if (this->hasAuxStorage(operation)) { - return false; + bool LinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const { + if (!auxiliaryRepeatedMultiplyMemory) { + auxiliaryRepeatedMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixColumnCount()); + return true; } - - auxiliaryRepeatedMultiplyStorage = std::make_unique<std::vector<ValueType>>(this->getMatrixColumnCount()); - return true; + return false; } template<typename ValueType> - bool LinearEquationSolver<ValueType>::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool LinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const { if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { - bool result = auxiliaryRepeatedMultiplyStorage != nullptr; - if (result) { - auxiliaryRepeatedMultiplyStorage.reset(); + if (auxiliaryRepeatedMultiplyMemory) { + auxiliaryRepeatedMultiplyMemory.reset(); + return true; } - return result; } return false; } template<typename ValueType> - bool LinearEquationSolver<ValueType>::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool LinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { - if (auxiliaryRepeatedMultiplyStorage != nullptr) { - result = auxiliaryRepeatedMultiplyStorage->size() != this->getMatrixColumnCount(); - auxiliaryRepeatedMultiplyStorage->resize(this->getMatrixColumnCount()); + if (auxiliaryRepeatedMultiplyMemory) { + result = auxiliaryRepeatedMultiplyMemory->size() != this->getMatrixColumnCount(); + auxiliaryRepeatedMultiplyMemory->resize(this->getMatrixColumnCount()); } } return result; } template<typename ValueType> - bool LinearEquationSolver<ValueType>::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool LinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const { if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { - return auxiliaryRepeatedMultiplyStorage != nullptr; + return static_cast<bool>(auxiliaryRepeatedMultiplyMemory); } return false; } diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h index 2e18ecd0c..aa243d593 100644 --- a/src/solver/LinearEquationSolver.h +++ b/src/solver/LinearEquationSolver.h @@ -69,36 +69,36 @@ namespace storm { // 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. + * 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 storage was allocated. + * @return True iff auxiliary memory was allocated. */ - virtual bool allocateAuxStorage(LinearEquationSolverOperation operation); + virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const; /*! - * Destroys previously allocated auxiliary storage for the provided operation. + * Destroys previously allocated auxiliary memory for the provided operation. * - * @return True iff auxiliary storage was deallocated. + * @return True iff auxiliary memory was deallocated. */ - virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation); + virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const; /*! - * If the matrix dimensions changed and auxiliary storage was allocated, this function needs to be called to - * update the auxiliary storage. + * 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 storage was reallocated. + * @return True iff the auxiliary memory was reallocated. */ - virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation); + virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const; /*! - * Checks whether the solver has allocated auxiliary storage for the provided operation. + * Checks whether the solver has allocated auxiliary memory for the provided operation. * - * @return True iff auxiliary storage was previously allocated (and not yet deallocated). + * @return True iff auxiliary memory was previously allocated (and not yet deallocated). */ - virtual bool hasAuxStorage(LinearEquationSolverOperation operation) const; + virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const; private: /*! @@ -111,8 +111,8 @@ namespace storm { */ virtual uint64_t getMatrixColumnCount() const = 0; - // Auxiliary storage for repeated matrix-vector multiplication. - mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRepeatedMultiplyStorage; + // 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 bb070c036..a92f19dd9 100644 --- a/src/solver/MinMaxLinearEquationSolver.cpp +++ b/src/solver/MinMaxLinearEquationSolver.cpp @@ -80,17 +80,17 @@ namespace storm { } template<typename ValueType> - bool MinMaxLinearEquationSolver<ValueType>::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool MinMaxLinearEquationSolver<ValueType>::allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const { return false; } template<typename ValueType> - bool MinMaxLinearEquationSolver<ValueType>::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool MinMaxLinearEquationSolver<ValueType>::deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const { return false; } template<typename ValueType> - bool MinMaxLinearEquationSolver<ValueType>::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + bool MinMaxLinearEquationSolver<ValueType>::hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const { return false; } diff --git a/src/solver/MinMaxLinearEquationSolver.h b/src/solver/MinMaxLinearEquationSolver.h index 565306756..0582269fb 100644 --- a/src/solver/MinMaxLinearEquationSolver.h +++ b/src/solver/MinMaxLinearEquationSolver.h @@ -126,21 +126,21 @@ namespace storm { * * @return True iff auxiliary storage was allocated. */ - virtual bool allocateAuxStorage(MinMaxLinearEquationSolverOperation operation); + virtual bool allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const; /*! * Destroys previously allocated auxiliary storage for the provided operation. * * @return True iff auxiliary storage was deallocated. */ - virtual bool deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation); + 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 hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const; + 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. diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index 0eb5d758e..ea39a075c 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -84,12 +84,12 @@ namespace storm { } template<typename ValueType> - NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingStorage(nullptr) { + 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(nullptr), A(nullptr), settings(settings), auxiliarySolvingStorage(nullptr) { + 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)); } @@ -107,9 +107,9 @@ namespace storm { template<typename ValueType> void NativeLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const { - bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); + bool allocatedAuxStorage = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations); if (allocatedAuxStorage) { - auxiliarySolvingStorage = std::make_unique<std::vector<ValueType>>(x.size()); + this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations); } if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::GaussSeidel) { @@ -124,11 +124,11 @@ namespace storm { A->performSuccessiveOverRelaxationStep(omega, x, b); // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision<ValueType>(*auxiliarySolvingStorage, x, 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 backup the contents of x. if (!converged) { - *auxiliarySolvingStorage = x; + *auxiliarySolvingMemory = x; } // Increase iteration count so we can abort if convergence is too slow. @@ -139,7 +139,7 @@ namespace storm { std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A->getJacobiDecomposition(); std::vector<ValueType>* currentX = &x; - std::vector<ValueType>* nextX = auxiliarySolvingStorage.get(); + std::vector<ValueType>* nextX = auxiliarySolvingMemory.get(); // Set up additional environment variables. uint_fast64_t iterationCount = 0; @@ -163,14 +163,14 @@ 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 == auxiliarySolvingStorage.get()) { + if (currentX == auxiliarySolvingMemory.get()) { std::swap(x, *currentX); } } - // If we allocated auxiliary storage, we need to dispose of it now. + // If we allocated auxiliary memory, we need to dispose of it now. if (allocatedAuxStorage) { - auxiliarySolvingStorage.reset(); + this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations); } } @@ -202,51 +202,51 @@ namespace storm { } template<typename ValueType> - bool NativeLinearEquationSolver<ValueType>::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool NativeLinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - if (!auxiliarySolvingStorage) { - auxiliarySolvingStorage = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount()); + if (!auxiliarySolvingMemory) { + auxiliarySolvingMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount()); result = true; } } - result |= LinearEquationSolver<ValueType>::allocateAuxStorage(operation); + result |= LinearEquationSolver<ValueType>::allocateAuxMemory(operation); return result; } template<typename ValueType> - bool NativeLinearEquationSolver<ValueType>::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool NativeLinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliarySolvingStorage) { + if (auxiliarySolvingMemory) { result = true; + auxiliarySolvingMemory.reset(); } - auxiliarySolvingStorage.reset(); } - result |= LinearEquationSolver<ValueType>::deallocateAuxStorage(operation); + result |= LinearEquationSolver<ValueType>::deallocateAuxMemory(operation); return result; } template<typename ValueType> - bool NativeLinearEquationSolver<ValueType>::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool NativeLinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliarySolvingStorage) { - result = auxiliarySolvingStorage->size() != this->getMatrixColumnCount(); - auxiliarySolvingStorage->resize(this->getMatrixRowCount()); + if (auxiliarySolvingMemory) { + result = auxiliarySolvingMemory->size() != this->getMatrixColumnCount(); + auxiliarySolvingMemory->resize(this->getMatrixRowCount()); } } - result |= LinearEquationSolver<ValueType>::reallocateAuxStorage(operation); + result |= LinearEquationSolver<ValueType>::reallocateAuxMemory(operation); return result; } template<typename ValueType> - bool NativeLinearEquationSolver<ValueType>::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool NativeLinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - result |= auxiliarySolvingStorage != nullptr; + result |= static_cast<bool>(auxiliarySolvingMemory); } - result |= LinearEquationSolver<ValueType>::hasAuxStorage(operation); + result |= LinearEquationSolver<ValueType>::hasAuxMemory(operation); return result; } diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index 19722debe..436e5cd36 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -55,10 +55,10 @@ namespace storm { NativeLinearEquationSolverSettings<ValueType>& getSettings(); NativeLinearEquationSolverSettings<ValueType> 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; + 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; @@ -75,8 +75,8 @@ namespace storm { // The settings used by the solver. NativeLinearEquationSolverSettings<ValueType> settings; - // Auxiliary storage for the equation solving methods. - mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingStorage; + // 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 2af6527e7..4e0609e80 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/solver/StandardMinMaxLinearEquationSolver.cpp @@ -106,7 +106,7 @@ namespace storm { // 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); + solver->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations); Status status = Status::InProgress; uint64_t iterations = 0; @@ -183,14 +183,13 @@ namespace storm { template<typename ValueType> 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); - bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::SolveEquations); - if (allocatedAuxStorage) { - auxiliarySolvingMultiplyStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); - auxiliarySolvingVectorStorage = std::make_unique<std::vector<ValueType>>(x.size()); + bool allocatedAuxMemory = !this->hasAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations); + if (allocatedAuxMemory) { + this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations); } std::vector<ValueType>* currentX = &x; - std::vector<ValueType>* newX = auxiliarySolvingVectorStorage.get(); + 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; @@ -198,10 +197,10 @@ namespace storm { Status status = Status::InProgress; while (status == Status::InProgress) { // Compute x' = A*x + b. - solver->multiply(*currentX, &b, *auxiliarySolvingMultiplyStorage); + solver->multiply(*currentX, &b, *auxiliarySolvingMultiplyMemory); // Reduce the vector x' by applying min/max for all non-deterministic choices. - storm::utility::vector::reduceVectorMinOrMax(dir, *auxiliarySolvingMultiplyStorage, *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())) { @@ -218,37 +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 == auxiliarySolvingVectorStorage.get()) { + if (currentX == auxiliarySolvingVectorMemory.get()) { std::swap(x, *currentX); } - // If we allocated auxiliary storage, we need to dispose of it now. - if (allocatedAuxStorage) { - auxiliarySolvingMultiplyStorage.reset(); - auxiliarySolvingVectorStorage.reset(); + // 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) const { - bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly); - if (allocatedAuxStorage) { - auxiliaryRepeatedMultiplyStorage = std::make_unique<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, b, *auxiliaryRepeatedMultiplyStorage); + 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, *auxiliaryRepeatedMultiplyStorage, x, this->A.getRowGroupIndices()); + storm::utility::vector::reduceVectorMinOrMax(dir, *auxiliaryRepeatedMultiplyMemory, x, this->A.getRowGroupIndices()); } - // If we allocated auxiliary storage, we need to dispose of it now. - if (allocatedAuxStorage) { - auxiliaryRepeatedMultiplyStorage.reset(); + // If we allocated auxiliary memory, we need to dispose of it now. + if (allocatedAuxMemory) { + this->deallocateAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly); } } @@ -286,47 +284,44 @@ namespace storm { } template<typename ValueType> - bool StandardMinMaxLinearEquationSolver<ValueType>::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool StandardMinMaxLinearEquationSolver<ValueType>::allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const { bool result = false; if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) { - if (!auxiliarySolvingMultiplyStorage) { + if (!auxiliarySolvingMultiplyMemory) { result = true; - auxiliarySolvingMultiplyStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); + auxiliarySolvingMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); } - if (!auxiliarySolvingVectorStorage) { + if (!auxiliarySolvingVectorMemory) { result = true; - auxiliarySolvingVectorStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowGroupCount()); + auxiliarySolvingVectorMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowGroupCount()); } } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) { - if (!auxiliarySolvingVectorStorage) { - result = true; - auxiliarySolvingVectorStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowGroupCount()); - } + // Nothing to do in this case. } else { STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method."); } } else { - if (!auxiliaryRepeatedMultiplyStorage) { + if (!auxiliaryRepeatedMultiplyMemory) { result = true; - auxiliaryRepeatedMultiplyStorage = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); + auxiliaryRepeatedMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowCount()); } } return result; } template<typename ValueType> - bool StandardMinMaxLinearEquationSolver<ValueType>::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool StandardMinMaxLinearEquationSolver<ValueType>::deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const { bool result = false; if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) { - if (auxiliarySolvingMultiplyStorage) { + if (auxiliarySolvingMultiplyMemory) { result = true; - auxiliarySolvingMultiplyStorage.reset(); + auxiliarySolvingMultiplyMemory.reset(); } - if (auxiliarySolvingVectorStorage) { + if (auxiliarySolvingVectorMemory) { result = true; - auxiliarySolvingVectorStorage.reset(); + auxiliarySolvingVectorMemory.reset(); } } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) { // Nothing to do in this case. @@ -334,23 +329,24 @@ namespace storm { STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method."); } } else { - if (auxiliaryRepeatedMultiplyStorage) { + if (auxiliaryRepeatedMultiplyMemory) { result = true; - auxiliaryRepeatedMultiplyStorage.reset(); + auxiliaryRepeatedMultiplyMemory.reset(); } } + return result; } template<typename ValueType> - bool StandardMinMaxLinearEquationSolver<ValueType>::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + bool StandardMinMaxLinearEquationSolver<ValueType>::hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const { if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) { - return auxiliarySolvingMultiplyStorage && auxiliarySolvingVectorStorage; + return auxiliarySolvingMultiplyMemory && auxiliarySolvingVectorMemory; } else { return false; } } else { - return static_cast<bool>(auxiliaryRepeatedMultiplyStorage); + return static_cast<bool>(auxiliaryRepeatedMultiplyMemory); } } diff --git a/src/solver/StandardMinMaxLinearEquationSolver.h b/src/solver/StandardMinMaxLinearEquationSolver.h index 523266c04..0c0b37f79 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/solver/StandardMinMaxLinearEquationSolver.h @@ -39,14 +39,14 @@ namespace storm { 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) const override; - virtual void multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) 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 allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; - virtual bool deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; - virtual bool hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const override; + 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) const; @@ -75,13 +75,12 @@ namespace storm { // the reference refers to localA. storm::storage::SparseMatrix<ValueType> const& A; - // Auxiliary storage for equation solving. - // Auxiliary storage for repeated matrix-vector multiplication. - mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingMultiplyStorage; - mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingVectorStorage; + // Auxiliary memory for equation solving. + mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingMultiplyMemory; + mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingVectorMemory; - // Auxiliary storage for repeated matrix-vector multiplication. - mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRepeatedMultiplyStorage; + // 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()); }