From be5fdeb636ebc049242e5d2c63446e8f4c653588 Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 28 Jun 2016 15:54:45 +0200 Subject: [PATCH 1/3] started working on internal auxiliary storage of solvers Former-commit-id: d895041c50a04955b365c0a192f8fd10adc44dfb --- .../csl/helper/SparseCtmcCslHelper.cpp | 9 +- .../helper/SparseMarkovAutomatonCslHelper.cpp | 7 +- .../prctl/helper/HybridMdpPrctlHelper.cpp | 6 +- .../prctl/helper/SparseDtmcPrctlHelper.cpp | 2 +- .../prctl/helper/SparseMdpPrctlHelper.cpp | 12 +- src/solver/EigenLinearEquationSolver.cpp | 31 +++- src/solver/EigenLinearEquationSolver.h | 12 +- .../EliminationLinearEquationSolver.cpp | 42 +++-- src/solver/EliminationLinearEquationSolver.h | 16 +- src/solver/GmmxxLinearEquationSolver.cpp | 121 +++++++++--- src/solver/GmmxxLinearEquationSolver.h | 26 ++- src/solver/LinearEquationSolver.cpp | 72 ++++++-- src/solver/LinearEquationSolver.h | 72 +++++++- src/solver/MinMaxLinearEquationSolver.cpp | 23 ++- src/solver/MinMaxLinearEquationSolver.h | 45 +++-- src/solver/NativeLinearEquationSolver.cpp | 153 ++++++++++------ src/solver/NativeLinearEquationSolver.h | 24 ++- .../StandardMinMaxLinearEquationSolver.cpp | 172 ++++++++++++------ .../StandardMinMaxLinearEquationSolver.h | 20 +- 19 files changed, 634 insertions(+), 231 deletions(-) diff --git a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index 4ccf61e64..186762a3d 100644 --- a/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -627,19 +627,19 @@ namespace storm { result = std::vector(values.size()); } } - std::vector multiplicationResult(result.size()); std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(uniformizedMatrix)); + solver->allocateAuxStorage(storm::solver::LinearEquationSolverOperation::MultiplyRepeatedly); if (!useMixedPoissonProbabilities && std::get<0>(foxGlynnResult) > 1) { // Perform the matrix-vector multiplications (without adding). - solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1, &multiplicationResult); + solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1); } else if (useMixedPoissonProbabilities) { std::function addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; }; // For the iterations below the left truncation point, we need to add and scale the result with the uniformization rate. for (uint_fast64_t index = 1; index < startingIteration; ++index) { - solver->repeatedMultiply(values, nullptr, 1, &multiplicationResult); + solver->repeatedMultiply(values, nullptr, 1); storm::utility::vector::applyPointwise(result, values, result, addAndScale); } } @@ -649,7 +649,7 @@ namespace storm { ValueType weight = 0; std::function addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; }; for (uint_fast64_t index = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) { - solver->repeatedMultiply(values, addVector, 1, &multiplicationResult); + solver->repeatedMultiply(values, addVector, 1); weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)]; storm::utility::vector::applyPointwise(result, values, result, addAndScale); @@ -658,7 +658,6 @@ namespace storm { return result; } - template storm::storage::SparseMatrix SparseCtmcCslHelper::computeProbabilityMatrix(storm::storage::SparseMatrix const& rateMatrix, std::vector const& exitRates) { // Turn the rates into probabilities by scaling each row with the exit rate of the state. diff --git a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index e56b21d66..b77a19825 100644 --- a/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -86,6 +86,7 @@ namespace storm { } std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(aProbabilistic); + solver->allocateAuxStorage(storm::solver::MinMaxLinearEquationSolverOperation::SolveEquations); // Perform the actual value iteration // * loop until the step bound has been reached @@ -94,15 +95,13 @@ namespace storm { // and 1_G being the characteristic vector for all goal states. // * perform one timed-step using v_MS := A_MSwG * v_MS + A_MStoPS * v_PS + (A * 1_G)|MS std::vector markovianNonGoalValuesSwap(markovianNonGoalValues); - std::vector multiplicationResultScratchMemory(aProbabilistic.getRowCount()); - std::vector aProbabilisticScratchMemory(probabilisticNonGoalValues.size()); for (uint_fast64_t currentStep = 0; currentStep < numberOfSteps; ++currentStep) { // Start by (re-)computing bProbabilistic = bProbabilisticFixed + aProbabilisticToMarkovian * vMarkovian. aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic); storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic); // Now perform the inner value iteration for probabilistic states. - solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory); + solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic); // (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic. aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian); @@ -115,7 +114,7 @@ namespace storm { // After the loop, perform one more step of the value iteration for PS states. aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic); storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic); - solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory); + solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic); } template diff --git a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp index 4863c1daf..2f3111a7f 100644 --- a/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp @@ -141,7 +141,7 @@ namespace storm { std::pair, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first)); - solver->multiply(dir, x, &explicitRepresentation.second, stepBound); + solver->repeatedMultiply(dir, x, &explicitRepresentation.second, stepBound); // Return a hybrid check result that stores the numerical values explicitly. return std::unique_ptr(new storm::modelchecker::HybridQuantitativeCheckResult(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.template toAdd(), maybeStates, odd, x)); @@ -166,7 +166,7 @@ namespace storm { // Perform the matrix-vector multiplication. std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitMatrix)); - solver->multiply(dir, x, nullptr, stepBound); + solver->repeatedMultiply(dir, x, nullptr, stepBound); // Return a hybrid check result that stores the numerical values explicitly. return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), odd, x)); @@ -195,7 +195,7 @@ namespace storm { // Perform the matrix-vector multiplication. std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first)); - solver->multiply(dir, x, &explicitRepresentation.second, stepBound); + solver->repeatedMultiply(dir, x, &explicitRepresentation.second, stepBound); // Return a hybrid check result that stores the numerical values explicitly. return std::unique_ptr(new HybridQuantitativeCheckResult(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero(), model.getReachableStates(), odd, x)); diff --git a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 916d24d53..79a85fdfe 100644 --- a/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -121,7 +121,7 @@ namespace storm { // Perform one single matrix-vector multiplication. std::unique_ptr> solver = linearEquationSolverFactory.create(transitionMatrix); - solver->repeatedMultiply(result); + solver->repeatedMultiply(result, nullptr, 1); return result; } diff --git a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index 8c09b96b8..483169789 100644 --- a/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -49,7 +49,7 @@ namespace storm { std::vector subresult(maybeStates.getNumberOfSetBits()); std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(std::move(submatrix)); - solver->multiply(dir, subresult, &b, stepBound); + solver->repeatedMultiply(dir, subresult, &b, stepBound); // Set the values of the resulting vector accordingly. storm::utility::vector::setVectorValues(result, maybeStates, subresult); @@ -67,7 +67,7 @@ namespace storm { storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one()); std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix); - solver->multiply(dir, result); + solver->repeatedMultiply(dir, result); return result; } @@ -211,7 +211,7 @@ namespace storm { std::vector result(rewardModel.getStateRewardVector()); std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix); - solver->multiply(dir, result, nullptr, stepCount); + solver->repeatedMultiply(dir, result, nullptr, stepCount); return result; } @@ -235,7 +235,7 @@ namespace storm { } std::unique_ptr> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix); - solver->multiply(dir, result, &totalRewardVector, stepBound); + solver->repeatedMultiply(dir, result, &totalRewardVector, stepBound); return result; } @@ -601,13 +601,13 @@ namespace storm { if (fixedTargetStates.get(state)) { builder.addNextValue(currentRow, newGoalState, conditionProbabilities[state]); if (!storm::utility::isZero(conditionProbabilities[state])) { - builder.addNextValue(currentRow, newFailState, 1 - conditionProbabilities[state]); + builder.addNextValue(currentRow, newFailState, storm::utility::one() - conditionProbabilities[state]); } ++currentRow; } else if (conditionStates.get(state)) { builder.addNextValue(currentRow, newGoalState, targetProbabilities[state]); if (!storm::utility::isZero(targetProbabilities[state])) { - builder.addNextValue(currentRow, newStopState, 1 - targetProbabilities[state]); + builder.addNextValue(currentRow, newStopState, storm::utility::one() - targetProbabilities[state]); } ++currentRow; } else { diff --git a/src/solver/EigenLinearEquationSolver.cpp b/src/solver/EigenLinearEquationSolver.cpp index c00fe97b0..25c16f1ae 100644 --- a/src/solver/EigenLinearEquationSolver.cpp +++ b/src/solver/EigenLinearEquationSolver.cpp @@ -109,12 +109,23 @@ namespace storm { template EigenLinearEquationSolver::EigenLinearEquationSolver(storm::storage::SparseMatrix&& A, EigenLinearEquationSolverSettings const& settings) : settings(settings) { + this->setMatrix(std::move(A)); + } + + template + void EigenLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix(A); + } + + template + void EigenLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + // Take ownership of the matrix so it is destroyed after we have translated it to Eigen's format. storm::storage::SparseMatrix localA(std::move(A)); - eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix(localA); + this->setMatrix(localA); } template - void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { // Map the input vectors to Eigen's format. auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); auto eigenB = Eigen::Matrix::Map(b.data(), b.size()); @@ -197,7 +208,7 @@ namespace storm { } template - void EigenLinearEquationSolver::multiply(std::vector& x, std::vector& result, std::vector const* b) const { + void EigenLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { // Typedef the map-type so we don't have to spell it out. typedef decltype(Eigen::Matrix::Map(b->data(), b->size())) MapType; @@ -234,10 +245,20 @@ namespace storm { return settings; } + template + uint64_t EigenLinearEquationSolver::getMatrixRowCount() const { + return eigenA->rows(); + } + + template + uint64_t EigenLinearEquationSolver::getMatrixColumnCount() const { + return eigenA->cols(); + } + // Specialization form storm::RationalNumber template<> - void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { // Map the input vectors to Eigen's format. auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); auto eigenB = Eigen::Matrix::Map(b.data(), b.size()); @@ -250,7 +271,7 @@ namespace storm { // Specialization form storm::RationalFunction template<> - void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void EigenLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { // Map the input vectors to Eigen's format. auto eigenX = Eigen::Matrix::Map(x.data(), x.size()); auto eigenB = Eigen::Matrix::Map(b.data(), b.size()); diff --git a/src/solver/EigenLinearEquationSolver.h b/src/solver/EigenLinearEquationSolver.h index 43e34427e..04c6c9718 100644 --- a/src/solver/EigenLinearEquationSolver.h +++ b/src/solver/EigenLinearEquationSolver.h @@ -61,13 +61,19 @@ namespace storm { EigenLinearEquationSolver(storm::storage::SparseMatrix const& A, EigenLinearEquationSolverSettings const& settings = EigenLinearEquationSolverSettings()); EigenLinearEquationSolver(storm::storage::SparseMatrix&& A, EigenLinearEquationSolverSettings const& settings = EigenLinearEquationSolverSettings()); - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual void solveEquations(std::vector& x, std::vector const& b) const override; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; EigenLinearEquationSolverSettings& getSettings(); EigenLinearEquationSolverSettings const& getSettings() const; - + private: + virtual uint64_t getMatrixRowCount() const override; + virtual uint64_t getMatrixColumnCount() const override; + // The (eigen) matrix associated with this equation solver. std::unique_ptr> eigenA; diff --git a/src/solver/EliminationLinearEquationSolver.cpp b/src/solver/EliminationLinearEquationSolver.cpp index 75340d6ba..982d2a54b 100644 --- a/src/solver/EliminationLinearEquationSolver.cpp +++ b/src/solver/EliminationLinearEquationSolver.cpp @@ -34,19 +34,29 @@ namespace storm { } template - EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A, EliminationLinearEquationSolverSettings const& settings) : localA(nullptr), A(A), settings(settings) { - // Intentionally left empty. + EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A, EliminationLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings) { + this->setMatrix(A); } template - EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix&& A, EliminationLinearEquationSolverSettings const& settings) : localA(std::make_unique>(std::move(A))), A(*localA), settings(settings) { - // Intentionally left empty. + EliminationLinearEquationSolver::EliminationLinearEquationSolver(storm::storage::SparseMatrix&& A, EliminationLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings) { + this->setMatrix(std::move(A)); } template - void EliminationLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { - STORM_LOG_WARN_COND(multiplyResult == nullptr, "Providing scratch memory will not improve the performance of this solver."); - + void EliminationLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + this->A = &A; + localA.reset(); + } + + template + void EliminationLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + localA = std::make_unique>(std::move(A)); + this->A = localA.get(); + } + + template + void EliminationLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { // FIXME: This solver will not work for all input systems. More concretely, the current implementation will // not work for systems that have a 0 on the diagonal. This is not a restriction of this technique in general // but arbitrary matrices require pivoting, which is not currently implemented. @@ -59,7 +69,7 @@ namespace storm { if (localA) { localA->convertToEquationSystem(); } else { - locallyConvertedMatrix = A; + locallyConvertedMatrix = *A; locallyConvertedMatrix.convertToEquationSystem(); } @@ -101,16 +111,16 @@ namespace storm { } template - void EliminationLinearEquationSolver::multiply(std::vector& x, std::vector& result, std::vector const* b) const { + void EliminationLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { if (&x != &result) { - A.multiplyWithVector(x, result); + A->multiplyWithVector(x, result); if (b != nullptr) { storm::utility::vector::addVectors(result, *b, result); } } else { // If the two vectors are aliases, we need to create a temporary. std::vector tmp(result.size()); - A.multiplyWithVector(x, tmp); + A->multiplyWithVector(x, tmp); if (b != nullptr) { storm::utility::vector::addVectors(tmp, *b, result); } @@ -127,6 +137,16 @@ namespace storm { return settings; } + template + uint64_t EliminationLinearEquationSolver::getMatrixRowCount() const { + return this->A->getRowCount(); + } + + template + uint64_t EliminationLinearEquationSolver::getMatrixColumnCount() const { + return this->A->getColumnCount(); + } + template std::unique_ptr> EliminationLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::make_unique>(matrix, settings); diff --git a/src/solver/EliminationLinearEquationSolver.h b/src/solver/EliminationLinearEquationSolver.h index 15a1e5c25..5963c5763 100644 --- a/src/solver/EliminationLinearEquationSolver.h +++ b/src/solver/EliminationLinearEquationSolver.h @@ -29,8 +29,11 @@ namespace storm { EliminationLinearEquationSolver(storm::storage::SparseMatrix const& A, EliminationLinearEquationSolverSettings const& settings = EliminationLinearEquationSolverSettings()); EliminationLinearEquationSolver(storm::storage::SparseMatrix&& A, EliminationLinearEquationSolverSettings const& settings = EliminationLinearEquationSolverSettings()); - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual void solveEquations(std::vector& x, std::vector const& b) const override; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; EliminationLinearEquationSolverSettings& getSettings(); EliminationLinearEquationSolverSettings const& getSettings() const; @@ -38,13 +41,16 @@ namespace storm { private: void initializeSettings(); + virtual uint64_t getMatrixRowCount() const override; + virtual uint64_t getMatrixColumnCount() const override; + // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted // when the solver is destructed. std::unique_ptr> localA; - // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix - // the reference refers to localA. - storm::storage::SparseMatrix const& A; + // A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix + // the pointer refers to localA. + storm::storage::SparseMatrix const* A; // The settings used by the solver. EliminationLinearEquationSolverSettings settings; diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index d29cda389..139a969dd 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -111,17 +111,31 @@ namespace storm { } template - GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(A), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A)), settings(settings) { - // Intentionally left empty. + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiStorage(nullptr) { + this->setMatrix(A); } template - GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings) : localA(std::make_unique>(std::move(A))), A(*localA), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(*localA)), settings(settings) { - // Intentionally left empty. + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiStorage(nullptr) { + this->setMatrix(std::move(A)); } template - void GmmxxLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void GmmxxLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + localA.reset(); + this->A = &A; + gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(A); + } + + template + void GmmxxLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + localA = std::make_unique>(std::move(A)); + this->A = localA.get(); + gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(*localA); + } + + template + void GmmxxLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { auto method = this->getSettings().getSolutionMethod(); auto preconditioner = this->getSettings().getPreconditioner(); STORM_LOG_INFO("Using method '" << method << "' with preconditioner '" << preconditioner << "' (max. " << this->getSettings().getMaximalNumberOfIterations() << " iterations)."); @@ -166,7 +180,7 @@ namespace storm { STORM_LOG_WARN("Iterative solver did not converge."); } } else if (method == GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi) { - uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b, multiplyResult); + uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(*A, x, b); // Check if the solver converged and issue a warning otherwise. if (iterations < this->getSettings().getMaximalNumberOfIterations()) { @@ -178,7 +192,7 @@ namespace storm { } template - void GmmxxLinearEquationSolver::multiply(std::vector& x, std::vector& result, std::vector const* b) const { + void GmmxxLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { if (b) { gmm::mult_add(*gmmxxMatrix, x, *b, result); } else { @@ -187,25 +201,20 @@ namespace storm { } template - uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { + bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliaryJacobiStorage = std::make_unique>(x.size()); + } + // Get a Jacobi decomposition of the matrix A. std::pair, std::vector> jacobiDecomposition = A.getJacobiDecomposition(); // Convert the LU matrix to gmm++'s format. std::unique_ptr> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.first)); - // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with. - bool multiplyResultProvided = true; - std::vector* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; std::vector* currentX = &x; - - // Target vector for precision calculation. - std::vector tmpX(x.size()); + std::vector* nextX = auxiliaryJacobiStorage.get(); // Set up additional environment variables. uint_fast64_t iterationCount = 0; @@ -213,9 +222,8 @@ namespace storm { while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) { // Compute D^-1 * (b - LU * x) and store result in nextX. - gmm::mult(*gmmLU, *currentX, tmpX); - gmm::add(b, gmm::scaled(tmpX, -storm::utility::one()), tmpX); - storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX); + gmm::mult_add(*gmmLU, gmm::scaled(*currentX, -storm::utility::one()), b, *nextX); + storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX); // Now check if the process already converged within our precision. converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion()); @@ -229,13 +237,13 @@ namespace storm { // If the last iteration did not write to the original x we have to swap the contents, because the // output has to be written to the input parameter x. - if (currentX == copyX) { + if (currentX == auxiliaryJacobiStorage.get()) { std::swap(x, *currentX); } - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; + // If we allocated auxiliary storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliaryJacobiStorage.reset(); } return iterationCount; @@ -251,6 +259,67 @@ namespace storm { return settings; } + template + bool GmmxxLinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi) { + if (!auxiliaryJacobiStorage) { + auxiliaryJacobiStorage = std::make_unique>(this->getMatrixRowCount()); + result = true; + } + } + } + result |= LinearEquationSolver::allocateAuxStorage(operation); + return result; + } + + template + bool GmmxxLinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliaryJacobiStorage) { + result = true; + } + auxiliaryJacobiStorage.reset(); + } + result |= LinearEquationSolver::deallocateAuxStorage(operation); + return result; + } + + template + bool GmmxxLinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliaryJacobiStorage != nullptr) { + result = auxiliaryJacobiStorage->size() != this->getMatrixColumnCount(); + auxiliaryJacobiStorage->resize(this->getMatrixRowCount()); + } + } + result |= LinearEquationSolver::reallocateAuxStorage(operation); + return result; + } + + template + bool GmmxxLinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + result |= auxiliaryJacobiStorage != nullptr; + } + result |= LinearEquationSolver::hasAuxStorage(operation); + return result; + } + + template + uint64_t GmmxxLinearEquationSolver::getMatrixRowCount() const { + return this->A->getRowCount(); + } + + template + uint64_t GmmxxLinearEquationSolver::getMatrixColumnCount() const { + return this->A->getColumnCount(); + } + template std::unique_ptr> GmmxxLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::make_unique>(matrix, settings); diff --git a/src/solver/GmmxxLinearEquationSolver.h b/src/solver/GmmxxLinearEquationSolver.h index 4425ec940..57386677d 100644 --- a/src/solver/GmmxxLinearEquationSolver.h +++ b/src/solver/GmmxxLinearEquationSolver.h @@ -89,12 +89,20 @@ namespace storm { GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings = GmmxxLinearEquationSolverSettings()); GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings = GmmxxLinearEquationSolverSettings()); - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual void solveEquations(std::vector& x, std::vector const& b) const override; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; GmmxxLinearEquationSolverSettings& getSettings(); GmmxxLinearEquationSolverSettings const& getSettings() const; + virtual bool allocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool hasAuxStorage(LinearEquationSolverOperation operation) const override; + private: /*! * Solves the linear equation system A*x = b given by the parameters using the Jacobi method. @@ -108,21 +116,27 @@ namespace storm { * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this * vector must be equal to the number of rows of A. */ - uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const; + uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const; + + virtual uint64_t getMatrixRowCount() const override; + virtual uint64_t getMatrixColumnCount() const override; // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted // when the solver is destructed. std::unique_ptr> localA; - // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix - // the reference refers to localA. - storm::storage::SparseMatrix const& A; + // A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix + // the pointer refers to localA. + storm::storage::SparseMatrix const* A; // The (gmm++) matrix associated with this equation solver. std::unique_ptr> gmmxxMatrix; // The settings used by the solver. GmmxxLinearEquationSolverSettings settings; + + // Auxiliary storage for the Jacobi method. + mutable std::unique_ptr> auxiliaryJacobiStorage; }; template diff --git a/src/solver/LinearEquationSolver.cpp b/src/solver/LinearEquationSolver.cpp index 84f5f84e5..feafdfde2 100644 --- a/src/solver/LinearEquationSolver.cpp +++ b/src/solver/LinearEquationSolver.cpp @@ -14,36 +14,80 @@ namespace storm { namespace solver { template - void LinearEquationSolver::repeatedMultiply(std::vector& x, std::vector const* b, uint_fast64_t n, std::vector* multiplyResult) const { + LinearEquationSolver::LinearEquationSolver() : auxiliaryRepeatedMultiplyStorage(nullptr) { + // Intentionally left empty. + } + + template + void LinearEquationSolver::repeatedMultiply(std::vector& x, std::vector const* b, uint_fast64_t n) const { + bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::MultiplyRepeatedly); + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage = std::make_unique>(x.size()); + } // Set up some temporary variables so that we can just swap pointers instead of copying the result after // each iteration. std::vector* currentX = &x; - - bool multiplyResultProvided = true; - std::vector* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; + std::vector* nextX = auxiliaryRepeatedMultiplyStorage.get(); // Now perform matrix-vector multiplication as long as we meet the bound. for (uint_fast64_t i = 0; i < n; ++i) { - this->multiply(*currentX, *nextX, b); + this->multiply(*currentX, b, *nextX); std::swap(nextX, currentX); } // If we performed an odd number of repetitions, we need to swap the contents of currentVector and x, // because the output is supposed to be stored in the input vector x. - if (currentX == copyX) { + if (currentX == auxiliaryRepeatedMultiplyStorage.get()) { std::swap(x, *currentX); } + + // If we allocated auxiliary storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage.reset(); + } + } + + template + bool LinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { + if (this->hasAuxStorage(operation)) { + return false; + } - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; + auxiliaryRepeatedMultiplyStorage = std::make_unique>(this->getMatrixColumnCount()); + return true; + } + + template + bool LinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { + bool result = auxiliaryRepeatedMultiplyStorage != nullptr; + if (result) { + auxiliaryRepeatedMultiplyStorage.reset(); + } + return result; + } + return false; + } + + template + bool LinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { + if (auxiliaryRepeatedMultiplyStorage != nullptr) { + result = auxiliaryRepeatedMultiplyStorage->size() != this->getMatrixColumnCount(); + auxiliaryRepeatedMultiplyStorage->resize(this->getMatrixColumnCount()); + } + } + return result; + } + + template + bool LinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { + return auxiliaryRepeatedMultiplyStorage != nullptr; } + return false; } template diff --git a/src/solver/LinearEquationSolver.h b/src/solver/LinearEquationSolver.h index 9632ddb2b..2e18ecd0c 100644 --- a/src/solver/LinearEquationSolver.h +++ b/src/solver/LinearEquationSolver.h @@ -11,6 +11,10 @@ namespace storm { namespace solver { + enum class LinearEquationSolverOperation { + SolveEquations, MultiplyRepeatedly + }; + /*! * An interface that represents an abstract linear equation solver. In addition to solving a system of linear * equations, the functionality to repeatedly multiply a matrix with a given vector is provided. @@ -18,10 +22,15 @@ namespace storm { template class LinearEquationSolver : public AbstractEquationSolver { public: + LinearEquationSolver(); + virtual ~LinearEquationSolver() { // Intentionally left empty. } + virtual void setMatrix(storm::storage::SparseMatrix const& A) = 0; + virtual void setMatrix(storm::storage::SparseMatrix&& A) = 0; + /*! * Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution. * The solution of the set of linear equations will be written to the vector x. Note that the matrix A has @@ -29,22 +38,20 @@ namespace storm { * * @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A. * @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A. - * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this - * vector must be equal to the number of rows of A. */ - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const = 0; + virtual void solveEquations(std::vector& x, std::vector const& b) const = 0; /*! * Performs on matrix-vector multiplication x' = A*x + b. * * @param x The input vector with which to multiply the matrix. Its length must be equal * to the number of columns of A. - * @param result The target vector into which to write the multiplication result. Its length must be equal - * to the number of rows of A. * @param b If non-null, this vector is added after the multiplication. If given, its length must be equal * to the number of rows of A. + * @param result The target vector into which to write the multiplication result. Its length must be equal + * to the number of rows of A. */ - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const = 0; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const = 0; /*! * Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After @@ -55,10 +62,57 @@ namespace storm { * to the number of columns of A. * @param b If non-null, this vector is added after each multiplication. If given, its length must be equal * to the number of rows of A. - * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this - * vector must be equal to the number of rows of A. + * @param n The number of times to perform the multiplication. + */ + void repeatedMultiply(std::vector& x, std::vector const* b, uint_fast64_t n) const; + + // Methods related to allocating/freeing auxiliary storage. + + /*! + * Allocates auxiliary storage that can be used to perform the provided operation. Repeated calls to the + * corresponding function can then be run without allocating/deallocating this storage repeatedly. + * Note: Since the allocated storage is fit to the currently selected options of the solver, they must not + * be changed any more after allocating the auxiliary storage until the storage is deallocated again. + * + * @return True iff auxiliary storage was allocated. + */ + virtual bool allocateAuxStorage(LinearEquationSolverOperation operation); + + /*! + * Destroys previously allocated auxiliary storage for the provided operation. + * + * @return True iff auxiliary storage was deallocated. */ - void repeatedMultiply(std::vector& x, std::vector const* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const; + virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation); + + /*! + * If the matrix dimensions changed and auxiliary storage was allocated, this function needs to be called to + * update the auxiliary storage. + * + * @return True iff the auxiliary storage was reallocated. + */ + virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation); + + /*! + * Checks whether the solver has allocated auxiliary storage for the provided operation. + * + * @return True iff auxiliary storage was previously allocated (and not yet deallocated). + */ + virtual bool hasAuxStorage(LinearEquationSolverOperation operation) const; + + private: + /*! + * Retrieves the row count of the matrix associated with this solver. + */ + virtual uint64_t getMatrixRowCount() const = 0; + + /*! + * Retrieves the column count of the matrix associated with this solver. + */ + virtual uint64_t getMatrixColumnCount() const = 0; + + // Auxiliary storage for repeated matrix-vector multiplication. + mutable std::unique_ptr> auxiliaryRepeatedMultiplyStorage; }; template diff --git a/src/solver/MinMaxLinearEquationSolver.cpp b/src/solver/MinMaxLinearEquationSolver.cpp index bd522834f..bb070c036 100644 --- a/src/solver/MinMaxLinearEquationSolver.cpp +++ b/src/solver/MinMaxLinearEquationSolver.cpp @@ -28,15 +28,15 @@ namespace storm { } template - void MinMaxLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void MinMaxLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { STORM_LOG_THROW(isSet(this->direction), storm::exceptions::IllegalFunctionCallException, "Optimization direction not set."); - solveEquations(convert(this->direction), x, b, multiplyResult, newX); + solveEquations(convert(this->direction), x, b); } template - void MinMaxLinearEquationSolver::multiply( std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const { + void MinMaxLinearEquationSolver::repeatedMultiply( std::vector& x, std::vector* b, uint_fast64_t n) const { STORM_LOG_THROW(isSet(this->direction), storm::exceptions::IllegalFunctionCallException, "Optimization direction not set."); - return multiply(convert(this->direction), x, b, n, multiplyResult); + return repeatedMultiply(convert(this->direction), x, b, n); } template @@ -79,6 +79,21 @@ namespace storm { return std::move(scheduler.get()); } + template + bool MinMaxLinearEquationSolver::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + return false; + } + + template + bool MinMaxLinearEquationSolver::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + return false; + } + + template + bool MinMaxLinearEquationSolver::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + return false; + } + template MinMaxLinearEquationSolverFactory::MinMaxLinearEquationSolverFactory(bool trackScheduler) : trackScheduler(trackScheduler) { // Intentionally left empty. diff --git a/src/solver/MinMaxLinearEquationSolver.h b/src/solver/MinMaxLinearEquationSolver.h index 358ef0be1..565306756 100644 --- a/src/solver/MinMaxLinearEquationSolver.h +++ b/src/solver/MinMaxLinearEquationSolver.h @@ -20,6 +20,10 @@ namespace storm { namespace solver { + enum class MinMaxLinearEquationSolverOperation { + SolveEquations, MultiplyRepeatedly + }; + /*! * A class representing the interface that all min-max linear equation solvers shall implement. */ @@ -40,20 +44,15 @@ namespace storm { * @param x The solution vector x. The initial values of x represent a guess of the real values to the * solver, but may be ignored. * @param b The vector to add after matrix-vector multiplication. - * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this - * vector must be equal to the number of rows of A. - * @param newX If non-null, this memory is used as a scratch memory. If given, the length of this - * vector must be equal to the length of the vector x (and thus to the number of columns of A). - * @return The solution vector x of the system of linear equations as the content of the parameter x. */ - virtual void solveEquations(OptimizationDirection d, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const = 0; + virtual void solveEquations(OptimizationDirection d, std::vector& x, std::vector const& b) const = 0; /*! * Behaves the same as the other variant of solveEquations, with the distinction that * instead of providing the optimization direction as an argument, the internally set optimization direction * is used. Note: this method can only be called after setting the optimization direction. */ - void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const; + void solveEquations(std::vector& x, std::vector const& b) const; /*! * Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes @@ -68,18 +67,16 @@ namespace storm { * i.e. after the method returns, this vector will contain the computed values. * @param b If not null, this vector is added after each multiplication. * @param n Specifies the number of iterations the matrix-vector multiplication is performed. - * @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this - * vector must be equal to the number of rows of A. * @return The result of the repeated matrix-vector multiplication as the content of the vector x. */ - virtual void multiply(OptimizationDirection d, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const = 0; + virtual void repeatedMultiply(OptimizationDirection d, std::vector& x, std::vector* b, uint_fast64_t n = 1) const = 0; /*! * Behaves the same as the other variant of multiply, with the * distinction that instead of providing the optimization direction as an argument, the internally set * optimization direction is used. Note: this method can only be called after setting the optimization direction. */ - virtual void multiply( std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const; + virtual void repeatedMultiply(std::vector& x, std::vector* b , uint_fast64_t n) const; /*! * Sets an optimization direction to use for calls to methods that do not explicitly provide one. @@ -119,6 +116,32 @@ namespace storm { */ std::unique_ptr getScheduler(); + // Methods related to allocating/freeing auxiliary storage. + + /*! + * Allocates auxiliary storage that can be used to perform the provided operation. Repeated calls to the + * corresponding function can then be run without allocating/deallocating this storage repeatedly. + * Note: Since the allocated storage is fit to the currently selected options of the solver, they must not + * be changed any more after allocating the auxiliary storage until the storage is deallocated again. + * + * @return True iff auxiliary storage was allocated. + */ + virtual bool allocateAuxStorage(MinMaxLinearEquationSolverOperation operation); + + /*! + * Destroys previously allocated auxiliary storage for the provided operation. + * + * @return True iff auxiliary storage was deallocated. + */ + virtual bool deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation); + + /*! + * Checks whether the solver has allocated auxiliary storage for the provided operation. + * + * @return True iff auxiliary storage was previously allocated (and not yet deallocated). + */ + virtual bool hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const; + protected: /// The optimization direction to use for calls to functions that do not provide it explicitly. Can also be unset. OptimizationDirectionSetting direction; diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index b24c9dde4..0eb5d758e 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -84,70 +84,62 @@ namespace storm { } template - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(A), settings(settings) { - // Intentionally left empty. + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingStorage(nullptr) { + this->setMatrix(A); } template - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings) : localA(std::make_unique>(std::move(A))), A(*localA), settings(settings) { - // Intentionally left empty. + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingStorage(nullptr) { + this->setMatrix(std::move(A)); } template - void NativeLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { + void NativeLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + localA.reset(); + this->A = &A; + } + + template + void NativeLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + localA = std::make_unique>(std::move(A)); + this->A = localA.get(); + } + + template + void NativeLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { + bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliarySolvingStorage = std::make_unique>(x.size()); + } + if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::GaussSeidel) { // Define the omega used for SOR. ValueType omega = this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR ? this->getSettings().getOmega() : storm::utility::one(); - // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with. - bool tmpXProvided = true; - std::vector* tmpX = multiplyResult; - if (multiplyResult == nullptr) { - tmpX = new std::vector(x); - tmpXProvided = false; - } else { - *tmpX = x; - } - // Set up additional environment variables. uint_fast64_t iterationCount = 0; bool converged = false; while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations()) { - A.performSuccessiveOverRelaxationStep(omega, x, b); + A->performSuccessiveOverRelaxationStep(omega, x, b); // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision(x, *tmpX, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)); + converged = storm::utility::vector::equalModuloPrecision(*auxiliarySolvingStorage, x, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)); - // If we did not yet converge, we need to copy the contents of x to *tmpX. + // If we did not yet converge, we need to backup the contents of x. if (!converged) { - *tmpX = x; + *auxiliarySolvingStorage = x; } // Increase iteration count so we can abort if convergence is too slow. ++iterationCount; } - - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!tmpXProvided) { - delete tmpX; - } } else { // Get a Jacobi decomposition of the matrix A. - std::pair, std::vector> jacobiDecomposition = A.getJacobiDecomposition(); + std::pair, std::vector> jacobiDecomposition = A->getJacobiDecomposition(); - // To avoid copying the contents of the vector in the loop, we create a temporary x to swap with. - bool multiplyResultProvided = true; - std::vector* nextX = multiplyResult; - if (nextX == nullptr) { - nextX = new std::vector(x.size()); - multiplyResultProvided = false; - } - std::vector const* copyX = nextX; std::vector* currentX = &x; - - // Target vector for multiplication. - std::vector tmpX(x.size()); + std::vector* nextX = auxiliarySolvingStorage.get(); // Set up additional environment variables. uint_fast64_t iterationCount = 0; @@ -155,15 +147,15 @@ namespace storm { while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) { // Compute D^-1 * (b - LU * x) and store result in nextX. - jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX); - storm::utility::vector::subtractVectors(b, tmpX, tmpX); - storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX); - - // Swap the two pointers as a preparation for the next iteration. - std::swap(nextX, currentX); + jacobiDecomposition.first.multiplyWithVector(*currentX, *nextX); + storm::utility::vector::subtractVectors(b, *nextX, *nextX); + storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX); // Now check if the process already converged within our precision. converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()); + + // Swap the two pointers as a preparation for the next iteration. + std::swap(nextX, currentX); // Increase iteration count so we can abort if convergence is too slow. ++iterationCount; @@ -171,28 +163,28 @@ namespace storm { // If the last iteration did not write to the original x we have to swap the contents, because the // output has to be written to the input parameter x. - if (currentX == copyX) { + if (currentX == auxiliarySolvingStorage.get()) { std::swap(x, *currentX); } - - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; - } + } + + // If we allocated auxiliary storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliarySolvingStorage.reset(); } } template - void NativeLinearEquationSolver::multiply(std::vector& x, std::vector& result, std::vector const* b) const { + void NativeLinearEquationSolver::multiply(std::vector& x, std::vector const* b, std::vector& result) const { if (&x != &result) { - A.multiplyWithVector(x, result); + A->multiplyWithVector(x, result); if (b != nullptr) { storm::utility::vector::addVectors(result, *b, result); } } else { // If the two vectors are aliases, we need to create a temporary. std::vector tmp(result.size()); - A.multiplyWithVector(x, tmp); + A->multiplyWithVector(x, tmp); if (b != nullptr) { storm::utility::vector::addVectors(tmp, *b, result); } @@ -209,6 +201,65 @@ namespace storm { return settings; } + template + bool NativeLinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (!auxiliarySolvingStorage) { + auxiliarySolvingStorage = std::make_unique>(this->getMatrixRowCount()); + result = true; + } + } + result |= LinearEquationSolver::allocateAuxStorage(operation); + return result; + } + + template + bool NativeLinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliarySolvingStorage) { + result = true; + } + auxiliarySolvingStorage.reset(); + } + result |= LinearEquationSolver::deallocateAuxStorage(operation); + return result; + } + + template + bool NativeLinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + if (auxiliarySolvingStorage) { + result = auxiliarySolvingStorage->size() != this->getMatrixColumnCount(); + auxiliarySolvingStorage->resize(this->getMatrixRowCount()); + } + } + result |= LinearEquationSolver::reallocateAuxStorage(operation); + return result; + } + + template + bool NativeLinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool result = false; + if (operation == LinearEquationSolverOperation::SolveEquations) { + result |= auxiliarySolvingStorage != nullptr; + } + result |= LinearEquationSolver::hasAuxStorage(operation); + return result; + } + + template + uint64_t NativeLinearEquationSolver::getMatrixRowCount() const { + return this->A->getRowCount(); + } + + template + uint64_t NativeLinearEquationSolver::getMatrixColumnCount() const { + return this->A->getColumnCount(); + } + template std::unique_ptr> NativeLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { return std::make_unique>(matrix, settings); diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index 235985df6..19722debe 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -46,23 +46,37 @@ namespace storm { NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings = NativeLinearEquationSolverSettings()); NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings = NativeLinearEquationSolverSettings()); - virtual void solveEquations(std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr) const override; - virtual void multiply(std::vector& x, std::vector& result, std::vector const* b = nullptr) const override; + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual void solveEquations(std::vector& x, std::vector const& b) const override; + virtual void multiply(std::vector& x, std::vector const* b, std::vector& result) const override; NativeLinearEquationSolverSettings& getSettings(); NativeLinearEquationSolverSettings const& getSettings() const; + virtual bool allocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation) override; + virtual bool hasAuxStorage(LinearEquationSolverOperation operation) const override; + private: + virtual uint64_t getMatrixRowCount() const override; + virtual uint64_t getMatrixColumnCount() const override; + // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted // when the solver is destructed. std::unique_ptr> localA; - // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix - // the reference refers to localA. - storm::storage::SparseMatrix const& A; + // A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix + // the pointer refers to localA. + storm::storage::SparseMatrix const* A; // The settings used by the solver. NativeLinearEquationSolverSettings settings; + + // Auxiliary storage for the equation solving methods. + mutable std::unique_ptr> auxiliarySolvingStorage; }; template diff --git a/src/solver/StandardMinMaxLinearEquationSolver.cpp b/src/solver/StandardMinMaxLinearEquationSolver.cpp index db507ef68..2af6527e7 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/solver/StandardMinMaxLinearEquationSolver.cpp @@ -84,45 +84,38 @@ namespace storm { } template - void StandardMinMaxLinearEquationSolver::solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void StandardMinMaxLinearEquationSolver::solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const { switch (this->getSettings().getSolutionMethod()) { - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: solveEquationsValueIteration(dir, x, b, multiplyResult, newX); break; - case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: solveEquationsPolicyIteration(dir, x, b, multiplyResult, newX); break; + case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration: solveEquationsValueIteration(dir, x, b); break; + case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: solveEquationsPolicyIteration(dir, x, b); break; } } template - void StandardMinMaxLinearEquationSolver::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { - - // Create scratch memory if none was provided. - bool multiplyResultMemoryProvided = multiplyResult != nullptr; - if (multiplyResult == nullptr) { - multiplyResult = new std::vector(this->A.getRowCount()); - } - + void StandardMinMaxLinearEquationSolver::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { // Create the initial scheduler. std::vector scheduler(this->A.getRowGroupCount()); // Create a vector for storing the right-hand side of the inner equation system. std::vector subB(this->A.getRowGroupCount()); - - // Create a vector that the inner equation solver can use as scratch memory. - std::vector deterministicMultiplyResult(this->A.getRowGroupCount()); + + // Resolve the nondeterminism according to the current scheduler. + storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(scheduler, true); + submatrix.convertToEquationSystem(); + storm::utility::vector::selectVectorValues(subB, scheduler, this->A.getRowGroupIndices(), b); + + // Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration. + auto solver = linearEquationSolverFactory->create(std::move(submatrix)); + solver->allocateAuxStorage(LinearEquationSolverOperation::SolveEquations); Status status = Status::InProgress; uint64_t iterations = 0; do { - // Resolve the nondeterminism according to the current scheduler. - storm::storage::SparseMatrix submatrix = this->A.selectRowsFromRowGroups(scheduler, true); - submatrix.convertToEquationSystem(); - storm::utility::vector::selectVectorValues(subB, scheduler, this->A.getRowGroupIndices(), b); - // Solve the equation system for the 'DTMC'. // FIXME: we need to remove the 0- and 1- states to make the solution unique. // HOWEVER: if we start with a valid scheduler, then we will never get an illegal one, because staying // within illegal MECs will never strictly improve the value. Is this true? - auto solver = linearEquationSolverFactory->create(std::move(submatrix)); - solver->solveEquations(x, subB, &deterministicMultiplyResult); + solver->solveEquations(x, subB); // Go through the multiplication result and see whether we can improve any of the choices. bool schedulerImproved = false; @@ -151,6 +144,12 @@ namespace storm { // If the scheduler did not improve, we are done. if (!schedulerImproved) { status = Status::Converged; + } else { + // Update the scheduler and the solver. + submatrix = this->A.selectRowsFromRowGroups(scheduler, true); + submatrix.convertToEquationSystem(); + storm::utility::vector::selectVectorValues(subB, scheduler, this->A.getRowGroupIndices(), b); + solver->setMatrix(std::move(submatrix)); } // Update environment variables. @@ -164,10 +163,6 @@ namespace storm { if (this->isTrackSchedulerSet()) { this->scheduler = std::make_unique(std::move(scheduler)); } - - if (!multiplyResultMemoryProvided) { - delete multiplyResult; - } } template @@ -186,22 +181,16 @@ namespace storm { } template - void StandardMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void StandardMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { std::unique_ptr> solver = linearEquationSolverFactory->create(A); - - // Create scratch memory if none was provided. - bool multiplyResultMemoryProvided = multiplyResult != nullptr; - if (multiplyResult == nullptr) { - multiplyResult = new std::vector(this->A.getRowCount()); - } - std::vector* currentX = &x; - bool xMemoryProvided = newX != nullptr; - if (newX == nullptr) { - newX = new std::vector(x.size()); + bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::SolveEquations); + if (allocatedAuxStorage) { + auxiliarySolvingMultiplyStorage = std::make_unique>(this->A.getRowCount()); + auxiliarySolvingVectorStorage = std::make_unique>(x.size()); } - // Keep track of which of the vectors for x is the auxiliary copy. - std::vector* copyX = newX; + std::vector* currentX = &x; + std::vector* newX = auxiliarySolvingVectorStorage.get(); // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations. uint64_t iterations = 0; @@ -209,10 +198,10 @@ namespace storm { Status status = Status::InProgress; while (status == Status::InProgress) { // Compute x' = A*x + b. - solver->multiply(*currentX, *multiplyResult, &b); + solver->multiply(*currentX, &b, *auxiliarySolvingMultiplyStorage); // Reduce the vector x' by applying min/max for all non-deterministic choices. - storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, *newX, this->A.getRowGroupIndices()); + storm::utility::vector::reduceVectorMinOrMax(dir, *auxiliarySolvingMultiplyStorage, *newX, this->A.getRowGroupIndices()); // Determine whether the method converged. if (storm::utility::vector::equalModuloPrecision(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) { @@ -229,39 +218,37 @@ namespace storm { // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result // is currently stored in currentX, but x is the output vector. - if (currentX == copyX) { + if (currentX == auxiliarySolvingVectorStorage.get()) { std::swap(x, *currentX); } - // Dispose of allocated scratch memory. - if (!xMemoryProvided) { - delete copyX; - } - if (!multiplyResultMemoryProvided) { - delete multiplyResult; + // If we allocated auxiliary storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliarySolvingMultiplyStorage.reset(); + auxiliarySolvingVectorStorage.reset(); } } template - void StandardMinMaxLinearEquationSolver::multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const { - std::unique_ptr> solver = linearEquationSolverFactory->create(A); - - // If scratch memory was not provided, we create it. - bool multiplyResultMemoryProvided = multiplyResult != nullptr; - if (!multiplyResult) { - multiplyResult = new std::vector(this->A.getRowCount()); + void StandardMinMaxLinearEquationSolver::multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const { + bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly); + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage = std::make_unique>(this->A.getRowCount()); } + std::unique_ptr> solver = linearEquationSolverFactory->create(A); + for (uint64_t i = 0; i < n; ++i) { - solver->multiply(x, *multiplyResult, b); + solver->multiply(x, b, *auxiliaryRepeatedMultiplyStorage); // Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost // element of the min/max operator stack. - storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices()); + storm::utility::vector::reduceVectorMinOrMax(dir, *auxiliaryRepeatedMultiplyStorage, x, this->A.getRowGroupIndices()); } - if (!multiplyResultMemoryProvided) { - delete multiplyResult; + // If we allocated auxiliary storage, we need to dispose of it now. + if (allocatedAuxStorage) { + auxiliaryRepeatedMultiplyStorage.reset(); } } @@ -298,6 +285,75 @@ namespace storm { return settings; } + template + bool StandardMinMaxLinearEquationSolver::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool result = false; + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { + if (!auxiliarySolvingMultiplyStorage) { + result = true; + auxiliarySolvingMultiplyStorage = std::make_unique>(this->A.getRowCount()); + } + if (!auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage = std::make_unique>(this->A.getRowGroupCount()); + } + } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration) { + if (!auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage = std::make_unique>(this->A.getRowGroupCount()); + } + } else { + STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method."); + } + } else { + if (!auxiliaryRepeatedMultiplyStorage) { + result = true; + auxiliaryRepeatedMultiplyStorage = std::make_unique>(this->A.getRowCount()); + } + } + return result; + } + + template + bool StandardMinMaxLinearEquationSolver::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool result = false; + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { + if (auxiliarySolvingMultiplyStorage) { + result = true; + auxiliarySolvingMultiplyStorage.reset(); + } + if (auxiliarySolvingVectorStorage) { + result = true; + auxiliarySolvingVectorStorage.reset(); + } + } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration) { + // Nothing to do in this case. + } else { + STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method."); + } + } else { + if (auxiliaryRepeatedMultiplyStorage) { + result = true; + auxiliaryRepeatedMultiplyStorage.reset(); + } + } + } + + template + bool StandardMinMaxLinearEquationSolver::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { + if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { + return auxiliarySolvingMultiplyStorage && auxiliarySolvingVectorStorage; + } else { + return false; + } + } else { + return static_cast(auxiliaryRepeatedMultiplyStorage); + } + } + template StandardMinMaxLinearEquationSolverFactory::StandardMinMaxLinearEquationSolverFactory(bool trackScheduler) : MinMaxLinearEquationSolverFactory(trackScheduler), linearEquationSolverFactory(nullptr) { // Intentionally left empty. diff --git a/src/solver/StandardMinMaxLinearEquationSolver.h b/src/solver/StandardMinMaxLinearEquationSolver.h index 36b2aa458..523266c04 100644 --- a/src/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/solver/StandardMinMaxLinearEquationSolver.h @@ -38,15 +38,19 @@ namespace storm { StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); - virtual void solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const override; - virtual void multiply(OptimizationDirection dir, std::vector& x, std::vector* b = nullptr, uint_fast64_t n = 1, std::vector* multiplyResult = nullptr) const override; + virtual void solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const override; + virtual void multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const override; StandardMinMaxLinearEquationSolverSettings const& getSettings() const; StandardMinMaxLinearEquationSolverSettings& getSettings(); + virtual bool allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; + virtual bool deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; + virtual bool hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const override; + private: - void solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const; - void solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const; + void solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + void solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const; @@ -70,6 +74,14 @@ namespace storm { // A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix // the reference refers to localA. storm::storage::SparseMatrix const& A; + + // Auxiliary storage for equation solving. + // Auxiliary storage for repeated matrix-vector multiplication. + mutable std::unique_ptr> auxiliarySolvingMultiplyStorage; + mutable std::unique_ptr> auxiliarySolvingVectorStorage; + + // Auxiliary storage for repeated matrix-vector multiplication. + mutable std::unique_ptr> auxiliaryRepeatedMultiplyStorage; }; template From 83c4b1647c7a9918926a82440e2eda551484813b Mon Sep 17 00:00:00 2001 From: dehnert Date: Tue, 28 Jun 2016 18:48:14 +0200 Subject: [PATCH 2/3] 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> 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> 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()); std::unique_ptr> 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 - GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiStorage(nullptr) { + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix const& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiMemory(nullptr) { this->setMatrix(A); } template - GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiStorage(nullptr) { + GmmxxLinearEquationSolver::GmmxxLinearEquationSolver(storm::storage::SparseMatrix&& A, GmmxxLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiMemory(nullptr) { this->setMatrix(std::move(A)); } @@ -202,9 +202,9 @@ namespace storm { template uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b) const { - bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); - if (allocatedAuxStorage) { - auxiliaryJacobiStorage = std::make_unique>(x.size()); + 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> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.first)); std::vector* currentX = &x; - std::vector* nextX = auxiliaryJacobiStorage.get(); + std::vector* 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 - bool GmmxxLinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool GmmxxLinearEquationSolver::allocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { if (this->getSettings().getSolutionMethod() == GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi) { - if (!auxiliaryJacobiStorage) { - auxiliaryJacobiStorage = std::make_unique>(this->getMatrixRowCount()); + if (!auxiliaryJacobiMemory) { + auxiliaryJacobiMemory = std::make_unique>(this->getMatrixRowCount()); result = true; } } } - result |= LinearEquationSolver::allocateAuxStorage(operation); + result |= LinearEquationSolver::allocateAuxMemory(operation); return result; } template - bool GmmxxLinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool GmmxxLinearEquationSolver::deallocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliaryJacobiStorage) { + if (auxiliaryJacobiMemory) { result = true; + auxiliaryJacobiMemory.reset(); } - auxiliaryJacobiStorage.reset(); } - result |= LinearEquationSolver::deallocateAuxStorage(operation); + result |= LinearEquationSolver::deallocateAuxMemory(operation); return result; } template - bool GmmxxLinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool GmmxxLinearEquationSolver::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::reallocateAuxStorage(operation); + result |= LinearEquationSolver::reallocateAuxMemory(operation); return result; } template - bool GmmxxLinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool GmmxxLinearEquationSolver::hasAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - result |= auxiliaryJacobiStorage != nullptr; + result |= static_cast(auxiliaryJacobiMemory); } - result |= LinearEquationSolver::hasAuxStorage(operation); + result |= LinearEquationSolver::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& getSettings(); GmmxxLinearEquationSolverSettings const& getSettings() const; - virtual bool allocateAuxStorage(LinearEquationSolverOperation operation) override; - virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation) override; - virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation) override; - virtual bool hasAuxStorage(LinearEquationSolverOperation operation) const override; + 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 settings; // Auxiliary storage for the Jacobi method. - mutable std::unique_ptr> auxiliaryJacobiStorage; + mutable std::unique_ptr> auxiliaryJacobiMemory; }; template 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 - LinearEquationSolver::LinearEquationSolver() : auxiliaryRepeatedMultiplyStorage(nullptr) { + LinearEquationSolver::LinearEquationSolver() : auxiliaryRepeatedMultiplyMemory(nullptr) { // Intentionally left empty. } template void LinearEquationSolver::repeatedMultiply(std::vector& x, std::vector const* b, uint_fast64_t n) const { - bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::MultiplyRepeatedly); - if (allocatedAuxStorage) { - auxiliaryRepeatedMultiplyStorage = std::make_unique>(x.size()); + 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* currentX = &x; - std::vector* nextX = auxiliaryRepeatedMultiplyStorage.get(); + std::vector* 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 - bool LinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { - if (this->hasAuxStorage(operation)) { - return false; + bool LinearEquationSolver::allocateAuxMemory(LinearEquationSolverOperation operation) const { + if (!auxiliaryRepeatedMultiplyMemory) { + auxiliaryRepeatedMultiplyMemory = std::make_unique>(this->getMatrixColumnCount()); + return true; } - - auxiliaryRepeatedMultiplyStorage = std::make_unique>(this->getMatrixColumnCount()); - return true; + return false; } template - bool LinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool LinearEquationSolver::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 - bool LinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool LinearEquationSolver::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 - bool LinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool LinearEquationSolver::hasAuxMemory(LinearEquationSolverOperation operation) const { if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) { - return auxiliaryRepeatedMultiplyStorage != nullptr; + return static_cast(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> auxiliaryRepeatedMultiplyStorage; + // Auxiliary memory for repeated matrix-vector multiplication. + mutable std::unique_ptr> auxiliaryRepeatedMultiplyMemory; }; template 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 - bool MinMaxLinearEquationSolver::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool MinMaxLinearEquationSolver::allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const { return false; } template - bool MinMaxLinearEquationSolver::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool MinMaxLinearEquationSolver::deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const { return false; } template - bool MinMaxLinearEquationSolver::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + bool MinMaxLinearEquationSolver::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 - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingStorage(nullptr) { + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingMemory(nullptr) { this->setMatrix(A); } template - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingStorage(nullptr) { + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix&& A, NativeLinearEquationSolverSettings const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingMemory(nullptr) { this->setMatrix(std::move(A)); } @@ -107,9 +107,9 @@ namespace storm { template void NativeLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { - bool allocatedAuxStorage = !this->hasAuxStorage(LinearEquationSolverOperation::SolveEquations); + bool allocatedAuxStorage = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations); if (allocatedAuxStorage) { - auxiliarySolvingStorage = std::make_unique>(x.size()); + this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations); } if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::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(*auxiliarySolvingStorage, x, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)); + converged = storm::utility::vector::equalModuloPrecision(*auxiliarySolvingMemory, x, static_cast(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, std::vector> jacobiDecomposition = A->getJacobiDecomposition(); std::vector* currentX = &x; - std::vector* nextX = auxiliarySolvingStorage.get(); + std::vector* 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 - bool NativeLinearEquationSolver::allocateAuxStorage(LinearEquationSolverOperation operation) { + bool NativeLinearEquationSolver::allocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - if (!auxiliarySolvingStorage) { - auxiliarySolvingStorage = std::make_unique>(this->getMatrixRowCount()); + if (!auxiliarySolvingMemory) { + auxiliarySolvingMemory = std::make_unique>(this->getMatrixRowCount()); result = true; } } - result |= LinearEquationSolver::allocateAuxStorage(operation); + result |= LinearEquationSolver::allocateAuxMemory(operation); return result; } template - bool NativeLinearEquationSolver::deallocateAuxStorage(LinearEquationSolverOperation operation) { + bool NativeLinearEquationSolver::deallocateAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - if (auxiliarySolvingStorage) { + if (auxiliarySolvingMemory) { result = true; + auxiliarySolvingMemory.reset(); } - auxiliarySolvingStorage.reset(); } - result |= LinearEquationSolver::deallocateAuxStorage(operation); + result |= LinearEquationSolver::deallocateAuxMemory(operation); return result; } template - bool NativeLinearEquationSolver::reallocateAuxStorage(LinearEquationSolverOperation operation) { + bool NativeLinearEquationSolver::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::reallocateAuxStorage(operation); + result |= LinearEquationSolver::reallocateAuxMemory(operation); return result; } template - bool NativeLinearEquationSolver::hasAuxStorage(LinearEquationSolverOperation operation) const { + bool NativeLinearEquationSolver::hasAuxMemory(LinearEquationSolverOperation operation) const { bool result = false; if (operation == LinearEquationSolverOperation::SolveEquations) { - result |= auxiliarySolvingStorage != nullptr; + result |= static_cast(auxiliarySolvingMemory); } - result |= LinearEquationSolver::hasAuxStorage(operation); + result |= LinearEquationSolver::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& getSettings(); NativeLinearEquationSolverSettings const& getSettings() const; - virtual bool allocateAuxStorage(LinearEquationSolverOperation operation) override; - virtual bool deallocateAuxStorage(LinearEquationSolverOperation operation) override; - virtual bool reallocateAuxStorage(LinearEquationSolverOperation operation) override; - virtual bool hasAuxStorage(LinearEquationSolverOperation operation) const override; + 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 settings; - // Auxiliary storage for the equation solving methods. - mutable std::unique_ptr> auxiliarySolvingStorage; + // Auxiliary memory for the equation solving methods. + mutable std::unique_ptr> auxiliarySolvingMemory; }; template 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 void StandardMinMaxLinearEquationSolver::solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const { std::unique_ptr> solver = linearEquationSolverFactory->create(A); - bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::SolveEquations); - if (allocatedAuxStorage) { - auxiliarySolvingMultiplyStorage = std::make_unique>(this->A.getRowCount()); - auxiliarySolvingVectorStorage = std::make_unique>(x.size()); + bool allocatedAuxMemory = !this->hasAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations); + if (allocatedAuxMemory) { + this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations); } std::vector* currentX = &x; - std::vector* newX = auxiliarySolvingVectorStorage.get(); + std::vector* 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(*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 - void StandardMinMaxLinearEquationSolver::multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const { - bool allocatedAuxStorage = !this->hasAuxStorage(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly); - if (allocatedAuxStorage) { - auxiliaryRepeatedMultiplyStorage = std::make_unique>(this->A.getRowCount()); + void StandardMinMaxLinearEquationSolver::repeatedMultiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const { + bool allocatedAuxMemory = !this->hasAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly); + if (allocatedAuxMemory) { + this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly); } std::unique_ptr> 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 - bool StandardMinMaxLinearEquationSolver::allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool StandardMinMaxLinearEquationSolver::allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const { bool result = false; if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { - if (!auxiliarySolvingMultiplyStorage) { + if (!auxiliarySolvingMultiplyMemory) { result = true; - auxiliarySolvingMultiplyStorage = std::make_unique>(this->A.getRowCount()); + auxiliarySolvingMultiplyMemory = std::make_unique>(this->A.getRowCount()); } - if (!auxiliarySolvingVectorStorage) { + if (!auxiliarySolvingVectorMemory) { result = true; - auxiliarySolvingVectorStorage = std::make_unique>(this->A.getRowGroupCount()); + auxiliarySolvingVectorMemory = std::make_unique>(this->A.getRowGroupCount()); } } else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration) { - if (!auxiliarySolvingVectorStorage) { - result = true; - auxiliarySolvingVectorStorage = std::make_unique>(this->A.getRowGroupCount()); - } + // 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>(this->A.getRowCount()); + auxiliaryRepeatedMultiplyMemory = std::make_unique>(this->A.getRowCount()); } } return result; } template - bool StandardMinMaxLinearEquationSolver::deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) { + bool StandardMinMaxLinearEquationSolver::deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const { bool result = false; if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::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::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 - bool StandardMinMaxLinearEquationSolver::hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const { + bool StandardMinMaxLinearEquationSolver::hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const { if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) { if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings::SolutionMethod::ValueIteration) { - return auxiliarySolvingMultiplyStorage && auxiliarySolvingVectorStorage; + return auxiliarySolvingMultiplyMemory && auxiliarySolvingVectorMemory; } else { return false; } } else { - return static_cast(auxiliaryRepeatedMultiplyStorage); + return static_cast(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&& A, std::unique_ptr>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings const& settings = StandardMinMaxLinearEquationSolverSettings()); virtual void solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const override; - virtual void multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const override; + virtual void repeatedMultiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const override; StandardMinMaxLinearEquationSolverSettings const& getSettings() const; StandardMinMaxLinearEquationSolverSettings& getSettings(); - virtual bool allocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; - virtual bool deallocateAuxStorage(MinMaxLinearEquationSolverOperation operation) override; - virtual bool hasAuxStorage(MinMaxLinearEquationSolverOperation operation) const override; + 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& x, std::vector const& b) const; @@ -75,13 +75,12 @@ namespace storm { // the reference refers to localA. storm::storage::SparseMatrix const& A; - // Auxiliary storage for equation solving. - // Auxiliary storage for repeated matrix-vector multiplication. - mutable std::unique_ptr> auxiliarySolvingMultiplyStorage; - mutable std::unique_ptr> auxiliarySolvingVectorStorage; + // Auxiliary memory for equation solving. + mutable std::unique_ptr> auxiliarySolvingMultiplyMemory; + mutable std::unique_ptr> auxiliarySolvingVectorMemory; - // Auxiliary storage for repeated matrix-vector multiplication. - mutable std::unique_ptr> auxiliaryRepeatedMultiplyStorage; + // Auxiliary memory for repeated matrix-vector multiplication. + mutable std::unique_ptr> auxiliaryRepeatedMultiplyMemory; }; template 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 - void TopologicalMinMaxLinearEquationSolver::solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult, std::vector* newX) const { + void TopologicalMinMaxLinearEquationSolver::solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const { #ifdef GPU_USE_FLOAT #define __FORCE_FLOAT_CALCULATION true @@ -49,7 +49,7 @@ namespace storm { std::vector new_x = storm::utility::vector::toValueType(x); std::vector const new_b = storm::utility::vector::toValueType(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 - void TopologicalMinMaxLinearEquationSolver::multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const { - - // If scratch memory was not provided, we need to create it. - bool multiplyResultMemoryProvided = true; - if (multiplyResult == nullptr) { - multiplyResult = new std::vector(this->A.getRowCount()); - multiplyResultMemoryProvided = false; - } + void TopologicalMinMaxLinearEquationSolver::repeatedMultiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const { + std::unique_ptr> multiplyResult = std::make_unique>(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 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 const& A, double precision = 1e-6, uint_fast64_t maximalNumberOfIterations = 20000, bool relative = true); - virtual void solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b, std::vector* multiplyResult = nullptr, std::vector* newX = nullptr) const override; + virtual void solveEquations(OptimizationDirection dir, std::vector& x, std::vector const& b) const override; - virtual void multiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n, std::vector* multiplyResult) const override; + virtual void repeatedMultiply(OptimizationDirection dir, std::vector& x, std::vector* b, uint_fast64_t n) const override; private: storm::storage::SparseMatrix 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(); 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().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().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().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().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().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(); 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().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().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().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().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().getPrecision()); } From f68120639399c1253f63872eae77b2434b5c24db Mon Sep 17 00:00:00 2001 From: TimQu Date: Wed, 29 Jun 2016 16:49:56 +0200 Subject: [PATCH 3/3] building markov automata from prism code Former-commit-id: 791c49c7cf4c2c0ab467afdd1ee3fcdc87433294 --- src/builder/ExplicitModelBuilder.cpp | 70 +++++++++++++++++-- src/builder/ExplicitModelBuilder.h | 19 ++++- src/generator/PrismNextStateGenerator.cpp | 4 +- src/models/sparse/MarkovAutomaton.cpp | 51 +++++++++++--- src/models/sparse/MarkovAutomaton.h | 41 ++++++++++- src/storage/BitVector.cpp | 10 +++ src/storage/BitVector.h | 11 +++ src/storage/SparseMatrix.cpp | 64 +++++++++++++++++ src/storage/SparseMatrix.h | 8 +++ src/utility/ModelInstantiator.h | 2 +- .../builder/ExplicitPrismModelBuilderTest.cpp | 26 +++++++ test/functional/builder/hybrid_states.ma | 19 +++++ test/functional/builder/simple.ma | 15 ++++ test/functional/builder/stream2.ma | 50 +++++++++++++ 14 files changed, 369 insertions(+), 21 deletions(-) create mode 100644 test/functional/builder/hybrid_states.ma create mode 100644 test/functional/builder/simple.ma create mode 100644 test/functional/builder/stream2.ma diff --git a/src/builder/ExplicitModelBuilder.cpp b/src/builder/ExplicitModelBuilder.cpp index f008321ee..aaa000000 100644 --- a/src/builder/ExplicitModelBuilder.cpp +++ b/src/builder/ExplicitModelBuilder.cpp @@ -5,6 +5,7 @@ #include "src/models/sparse/Dtmc.h" #include "src/models/sparse/Ctmc.h" #include "src/models/sparse/Mdp.h" +#include "src/models/sparse/MarkovAutomaton.h" #include "src/models/sparse/StandardRewardModel.h" #include "src/storage/expressions/ExpressionManager.h" @@ -133,6 +134,9 @@ namespace storm { case storm::generator::ModelType::MDP: result = std::shared_ptr>(new storm::models::sparse::Mdp(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), std::move(modelComponents.rewardModels), std::move(modelComponents.choiceLabeling))); break; + case storm::generator::ModelType::MA: + result = std::shared_ptr>(new storm::models::sparse::MarkovAutomaton(std::move(modelComponents.transitionMatrix), std::move(modelComponents.stateLabeling), *std::move(modelComponents.markovianStates), std::move(modelComponents.rewardModels), std::move(modelComponents.choiceLabeling))); + break; default: STORM_LOG_THROW(false, storm::exceptions::WrongFormatException, "Error while creating model: cannot handle this model type."); break; @@ -165,12 +169,17 @@ namespace storm { } template - boost::optional>> ExplicitModelBuilder::buildMatrices(storm::storage::SparseMatrixBuilder& transitionMatrixBuilder, std::vector>& rewardModelBuilders) { + void ExplicitModelBuilder::buildMatrices(storm::storage::SparseMatrixBuilder& transitionMatrixBuilder, std::vector>& rewardModelBuilders, boost::optional>>& choiceLabels , boost::optional& markovianChoices) { // Create choice labels, if requested, - boost::optional>> choiceLabels; if (generator->getOptions().isBuildChoiceLabelsSet()) { choiceLabels = std::vector>(); } + + // Create markovian states bit vector, if required + if (generator->getModelType() == storm::generator::ModelType::MA) { + // The BitVector will be resized when the correct size is known + markovianChoices = storm::storage::BitVector(64, false); + } // Create a callback for the next-state generator to enable it to request the index of states. std::function stateToIdCallback = std::bind(&ExplicitModelBuilder::getOrAddStateIndex, this, std::placeholders::_1); @@ -219,6 +228,12 @@ namespace storm { // Insert empty choice labeling for added self-loop transitions. choiceLabels.get().push_back(boost::container::flat_set()); } + + if (generator->getModelType() == storm::generator::ModelType::MA) { + markovianChoices->enlargeLiberally(currentRow, false); + markovianChoices->set(currentRow); + } + if (!generator->isDeterministicModel()) { transitionMatrixBuilder.newRowGroup(currentRow); } @@ -262,6 +277,12 @@ namespace storm { choiceLabels.get().push_back(choice.getChoiceLabels()); } + // If we keep track of the Markovian choices, store whether the current one is Markovian. + if( markovianChoices && choice.isMarkovian() ) { + markovianChoices->enlargeLiberally(currentRow, false); + markovianChoices->set(currentRow); + } + // Add the probabilistic behavior to the matrix. for (auto const& stateProbabilityPair : choice) { transitionMatrixBuilder.addNextValue(currentRow, stateProbabilityPair.first, stateProbabilityPair.second); @@ -280,6 +301,11 @@ namespace storm { ++currentRowGroup; } } + + if (markovianChoices) { + // We now know the correct size + markovianChoices->resize(currentRow, false); + } // If the exploration order was not breadth-first, we need to fix the entries in the matrix according to // (reversed) mapping of row groups to indices. @@ -304,8 +330,6 @@ namespace storm { // Fix (c). this->stateStorage.stateToId.remap([&remapping] (StateType const& state) { return remapping[state]; } ); } - - return choiceLabels; } template @@ -323,7 +347,8 @@ namespace storm { rewardModelBuilders.emplace_back(generator->getRewardModelInformation(i)); } - modelComponents.choiceLabeling = buildMatrices(transitionMatrixBuilder, rewardModelBuilders); + boost::optional markovianChoices; + buildMatrices(transitionMatrixBuilder, rewardModelBuilders, modelComponents.choiceLabeling, markovianChoices); modelComponents.transitionMatrix = transitionMatrixBuilder.build(); // Now finalize all reward models. @@ -342,9 +367,44 @@ namespace storm { } } + if (generator->getModelType() == storm::generator::ModelType::MA) { + STORM_LOG_ASSERT(markovianChoices, "No information regarding markovian choices available."); + buildMarkovianStates(modelComponents, *markovianChoices); + } + return modelComponents; } + template + void ExplicitModelBuilder::buildMarkovianStates(ModelComponents& modelComponents, storm::storage::BitVector const& markovianChoices) const { + modelComponents.markovianStates = storm::storage::BitVector(modelComponents.transitionMatrix.getRowGroupCount(), false); + // Check for each state whether it contains a markovian choice. + for (uint_fast64_t state = 0; state < modelComponents.transitionMatrix.getRowGroupCount(); ++state ) { + uint_fast64_t firstChoice = modelComponents.transitionMatrix.getRowGroupIndices()[state]; + uint_fast64_t markovianChoice = markovianChoices.getNextSetIndex(firstChoice); + if (markovianChoice < modelComponents.transitionMatrix.getRowGroupIndices()[state+1]) { + // Found a markovian choice. Assert that there is not a second one. + STORM_LOG_THROW(markovianChoices.getNextSetIndex(markovianChoice+1) >= modelComponents.transitionMatrix.getRowGroupIndices()[state+1], storm::exceptions::WrongFormatException, "Multiple Markovian choices defined for some state"); + modelComponents.markovianStates->set(state); + // Swap the first choice and the found markovian choice (if they are not equal) + if (firstChoice != markovianChoice) { + modelComponents.transitionMatrix.swapRows(firstChoice, markovianChoice); + for (auto& rewardModel : modelComponents.rewardModels) { + if (rewardModel.second.hasStateActionRewards()) { + std::swap(rewardModel.second.getStateActionRewardVector()[firstChoice], rewardModel.second.getStateActionRewardVector()[markovianChoice]); + } + if (rewardModel.second.hasTransitionRewards()) { + rewardModel.second.getTransitionRewardMatrix().swapRows(firstChoice, markovianChoice); + } + } + if (modelComponents.choiceLabeling) { + std::swap(modelComponents.choiceLabeling.get()[firstChoice], modelComponents.choiceLabeling.get()[markovianChoice]); + } + } + } + } + } + template storm::models::sparse::StateLabeling ExplicitModelBuilder::buildStateLabeling() { return generator->label(stateStorage.stateToId, stateStorage.initialStateIndices, stateStorage.deadlockStateIndices); diff --git a/src/builder/ExplicitModelBuilder.h b/src/builder/ExplicitModelBuilder.h index 296c15b7d..bb404f088 100644 --- a/src/builder/ExplicitModelBuilder.h +++ b/src/builder/ExplicitModelBuilder.h @@ -63,6 +63,9 @@ namespace storm { // A vector that stores a labeling for each choice. boost::optional>> choiceLabeling; + + // A vector that stores which states are markovian. + boost::optional markovianStates; }; struct Options { @@ -138,10 +141,10 @@ namespace storm { * * @param transitionMatrixBuilder The builder of the transition matrix. * @param rewardModelBuilders The builders for the selected reward models. - * @return A tuple containing a vector with all rows at which the nondeterministic choices of each state begin - * and a vector containing the labels associated with each choice. + * @param choiceLabels is set to a vector containing the labels associated with each choice (is only set if choice labels are requested). + * @param markovianChoices is set to a bit vector storing whether a choice is markovian (is only set if the model type requires this information). */ - boost::optional>> buildMatrices(storm::storage::SparseMatrixBuilder& transitionMatrixBuilder, std::vector>& rewardModelBuilders); + void buildMatrices(storm::storage::SparseMatrixBuilder& transitionMatrixBuilder, std::vector>& rewardModelBuilders, boost::optional>>& choiceLabels , boost::optional& markovianChoices); /*! * Explores the state space of the given program and returns the components of the model as a result. @@ -150,6 +153,16 @@ namespace storm { */ ModelComponents buildModelComponents(); + /*! + * Set the markovian states of the given modelComponents, + * makes sure that each state has at most one markovian choice, + * and makes this choice the first one of the corresponding state + * + * @param modelComponents The components of the model build so far + * @markovianChoices bit vector storing whether a choice is markovian + */ + void buildMarkovianStates(ModelComponents& modelComponents, storm::storage::BitVector const& markovianChoices) const; + /*! * Builds the state labeling for the given program. * diff --git a/src/generator/PrismNextStateGenerator.cpp b/src/generator/PrismNextStateGenerator.cpp index 1173fc5d2..552cd30a9 100644 --- a/src/generator/PrismNextStateGenerator.cpp +++ b/src/generator/PrismNextStateGenerator.cpp @@ -366,7 +366,7 @@ namespace storm { continue; } - result.push_back(Choice(command.getActionIndex())); + result.push_back(Choice(command.getActionIndex(), command.isMarkovian())); Choice& choice = result.back(); // Remember the command labels only if we were asked to. @@ -564,4 +564,4 @@ namespace storm { template class PrismNextStateGenerator; template class PrismNextStateGenerator; } -} \ No newline at end of file +} diff --git a/src/models/sparse/MarkovAutomaton.cpp b/src/models/sparse/MarkovAutomaton.cpp index 81f3bc5c9..f822de3c1 100644 --- a/src/models/sparse/MarkovAutomaton.cpp +++ b/src/models/sparse/MarkovAutomaton.cpp @@ -1,12 +1,14 @@ #include "src/models/sparse/MarkovAutomaton.h" -#include "src/models/sparse/StandardRewardModel.h" -#include "src/exceptions/InvalidArgumentException.h" -#include "src/utility/constants.h" + #include "src/adapters/CarlAdapter.h" -#include "src/storage/FlexibleSparseMatrix.h" -#include "src/models/sparse/Dtmc.h" +#include "src/models/sparse/StandardRewardModel.h" #include "src/solver/stateelimination/StateEliminator.h" +#include "src/storage/FlexibleSparseMatrix.h" +#include "src/utility/constants.h" #include "src/utility/vector.h" +#include "src/utility/macros.h" + +#include "src/exceptions/InvalidArgumentException.h" namespace storm { namespace models { @@ -33,6 +35,30 @@ namespace storm { : NondeterministicModel(storm::models::ModelType::MarkovAutomaton, std::move(transitionMatrix), std::move(stateLabeling), std::move(rewardModels), std::move(optionalChoiceLabeling)), markovianStates(markovianStates), exitRates(std::move(exitRates)), closed(this->checkIsClosed()) { this->turnRatesToProbabilities(); } + + + template + MarkovAutomaton::MarkovAutomaton(storm::storage::SparseMatrix const& rateMatrix, + storm::models::sparse::StateLabeling const& stateLabeling, + storm::storage::BitVector const& markovianStates, + std::unordered_map const& rewardModels, + boost::optional> const& optionalChoiceLabeling) + : NondeterministicModel(storm::models::ModelType::MarkovAutomaton, rateMatrix, stateLabeling, rewardModels, optionalChoiceLabeling), markovianStates(markovianStates) { + turnRatesToProbabilities(); + this->closed = checkIsClosed(); + } + + template + MarkovAutomaton::MarkovAutomaton(storm::storage::SparseMatrix&& rateMatrix, + storm::models::sparse::StateLabeling&& stateLabeling, + storm::storage::BitVector&& markovianStates, + std::unordered_map&& rewardModels, + boost::optional>&& optionalChoiceLabeling) + : NondeterministicModel(storm::models::ModelType::MarkovAutomaton, rateMatrix, stateLabeling, rewardModels, optionalChoiceLabeling), markovianStates(markovianStates) { + turnRatesToProbabilities(); + this->closed = checkIsClosed(); + } + template MarkovAutomaton::MarkovAutomaton(storm::storage::SparseMatrix&& transitionMatrix, @@ -230,9 +256,18 @@ namespace storm { template void MarkovAutomaton::turnRatesToProbabilities() { - for (auto state : this->markovianStates) { - for (auto& transition : this->getTransitionMatrix().getRow(this->getTransitionMatrix().getRowGroupIndices()[state])) { - transition.setValue(transition.getValue() / this->exitRates[state]); + this->exitRates.resize(this->getNumberOfStates()); + for (uint_fast64_t state = 0; state< this->getNumberOfStates(); ++state) { + uint_fast64_t row = this->getTransitionMatrix().getRowGroupIndices()[state]; + if(this->markovianStates.get(state)) { + this->exitRates[state] = this->getTransitionMatrix().getRowSum(row); + for (auto& transition : this->getTransitionMatrix().getRow(row)) { + transition.setValue(transition.getValue() / this->exitRates[state]); + } + ++row; + } + for(; row < this->getTransitionMatrix().getRowGroupIndices()[state+1]; ++row) { + STORM_LOG_THROW(storm::utility::isOne(this->getTransitionMatrix().getRowSum(row)), storm::exceptions::InvalidArgumentException, "Transitions of rateMatrix do not sum up to one for some non-Markovian choice."); } } } diff --git a/src/models/sparse/MarkovAutomaton.h b/src/models/sparse/MarkovAutomaton.h index 9864f4f85..f78dc073c 100644 --- a/src/models/sparse/MarkovAutomaton.h +++ b/src/models/sparse/MarkovAutomaton.h @@ -31,7 +31,43 @@ namespace storm { std::vector const& exitRates, std::unordered_map const& rewardModels = std::unordered_map(), boost::optional> const& optionalChoiceLabeling = boost::optional>()); + + /*! + * Constructs a model from the given data. + * + * For hybrid states (i.e., states with Markovian and probabilistic transitions), it is assumed that the first + * choice corresponds to the markovian transitions. + * + * @param rateMatrix The matrix representing the transitions in the model in terms of rates. + * @param stateLabeling The labeling of the states. + * @param markovianStates A bit vector indicating the Markovian states of the automaton (i.e., states with at least one markovian transition). + * @param rewardModels A mapping of reward model names to reward models. + * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state. + */ + MarkovAutomaton(storm::storage::SparseMatrix const& rateMatrix, + storm::models::sparse::StateLabeling const& stateLabeling, + storm::storage::BitVector const& markovianStates, + std::unordered_map const& rewardModels = std::unordered_map(), + boost::optional> const& optionalChoiceLabeling = boost::optional>()); + /*! + * Constructs a model from the given data. + * + * For hybrid states (i.e., states with Markovian and probabilistic transitions), it is assumed that the first + * choice corresponds to the markovian transitions. + * + * @param rateMatrix The matrix representing the transitions in the model in terms of rates. + * @param stateLabeling The labeling of the states. + * @param markovianStates A bit vector indicating the Markovian states of the automaton (i.e., states with at least one markovian transition). + * @param rewardModels A mapping of reward model names to reward models. + * @param optionalChoiceLabeling A vector that represents the labels associated with the choices of each state. + */ + MarkovAutomaton(storm::storage::SparseMatrix&& rateMatrix, + storm::models::sparse::StateLabeling&& stateLabeling, + storm::storage::BitVector&& markovianStates, + std::unordered_map&& rewardModels = std::unordered_map(), + boost::optional>&& optionalChoiceLabeling = boost::optional>()); + /*! * Constructs a model by moving the given data. * @@ -48,7 +84,7 @@ namespace storm { std::vector const& exitRates, std::unordered_map&& rewardModels = std::unordered_map(), boost::optional>&& optionalChoiceLabeling = boost::optional>()); - + /*! * Constructs a model by moving the given data. * @@ -160,7 +196,8 @@ namespace storm { /*! * Under the assumption that the Markovian choices of this Markov automaton are expressed in terms of * rates in the transition matrix, this procedure turns the rates into the corresponding probabilities by - * dividing each entry by the exit rate of the state. + * dividing each entry by the sum of the rates for that choice. + * Also sets the exitRates accordingly and throws an exception if the values for a non-markovian choice do not sum up to one. */ void turnRatesToProbabilities(); diff --git a/src/storage/BitVector.cpp b/src/storage/BitVector.cpp index a60b8c10c..ccbcce134 100644 --- a/src/storage/BitVector.cpp +++ b/src/storage/BitVector.cpp @@ -252,6 +252,16 @@ namespace storm { truncateLastBucket(); } } + + void BitVector::enlargeLiberally(uint_fast64_t minimumLength, bool init) { + if(minimumLength > this->size()) { + uint_fast64_t newLength = this->bucketCount() << 6; + while(newLength < minimumLength) { + newLength = newLength << 1; + } + resize(newLength, init); + } + } BitVector BitVector::operator&(BitVector const& other) const { STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match."); diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h index 62ae010be..6ae966379 100644 --- a/src/storage/BitVector.h +++ b/src/storage/BitVector.h @@ -234,6 +234,17 @@ namespace storm { */ void resize(uint_fast64_t newLength, bool init = false); + /*! + * Enlarges the bit vector such that it holds at least the given number of bits (but possibly more). + * This can be used to diminish reallocations when the final size of the bit vector is not known yet. + * The bit vector does not become smaller. + * New bits are initialized to the given value. + * + * @param minimumLength The minimum number of bits that the bit vector should hold. + * @param init The truth value to which to initialize newly created bits. + */ + void enlargeLiberally(uint_fast64_t minimumLength, bool init = false); + /*! * Performs a logical "and" with the given bit vector. In case the sizes of the bit vectors do not match, * only the matching portion is considered and the overlapping bits are set to 0. diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index e4e81e4c5..95f65ce99 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -658,6 +658,70 @@ namespace storm { return bv; } + template + void SparseMatrix::swapRows(index_type const& row1, index_type const& row2) { + if(row1==row2) { + return; + } + + // Get the index of the row that has more / less entries than the other + index_type largerRow = getRow(row1).getNumberOfEntries() > getRow(row2).getNumberOfEntries() ? row1 : row2; + index_type smallerRow = largerRow == row1 ? row2 : row1; + index_type rowSizeDifference = getRow(largerRow).getNumberOfEntries() - getRow(smallerRow).getNumberOfEntries(); + // Save contents of larger row + std::vector> largerRowContents(getRow(largerRow).begin(), getRow(largerRow).end()); + + if(largerRow < smallerRow) { + auto writeIt = getRows(largerRow, smallerRow+1).begin(); + // write smaller row in its new position + for(auto& smallerRowEntry : getRow(smallerRow)) { + *writeIt = std::move(smallerRowEntry); + ++writeIt; + } + if(!storm::utility::isZero(rowSizeDifference)) { + // write the intermediate rows into their correct position + for(auto& intermediateRowEntry : getRows(largerRow+1, smallerRow)) { + *writeIt = std::move(intermediateRowEntry); + ++writeIt; + } + } + // write the larger row + for(auto& largerRowEntry : largerRowContents) { + *writeIt = std::move(largerRowEntry); + ++writeIt; + } + STORM_LOG_ASSERT(writeIt == getRow(smallerRow).end(), "Unexpected position of write iterator"); + //Update row indications + for(index_type row = largerRow +1; row <= smallerRow; ++row) { + rowIndications[row] -= rowSizeDifference; + } + } else { + auto writeIt = getRows(smallerRow, largerRow+1).end() -1; + // write smaller row in its new position + for(auto smallerRowEntryIt = getRow(smallerRow).end() -1; smallerRowEntryIt != getRow(smallerRow).begin()-1; --smallerRowEntryIt) { + *writeIt = std::move(*smallerRowEntryIt); + --writeIt; + } + if(!storm::utility::isZero(rowSizeDifference)) { + // write the intermediate rows into their correct position + for(auto intermediateRowEntryIt = getRows(smallerRow+1, largerRow).end() -1; intermediateRowEntryIt != getRows(smallerRow+1, largerRow).begin()-1; --intermediateRowEntryIt) { + *writeIt = std::move(*intermediateRowEntryIt); + --writeIt; + } + } + // write the larger row + for(auto largerRowEntryIt = largerRowContents.rbegin(); largerRowEntryIt != largerRowContents.rend(); ++largerRowEntryIt) { + *writeIt = std::move(*largerRowEntryIt); + --writeIt; + } + STORM_LOG_ASSERT(writeIt == getRow(smallerRow).begin()-1, "Unexpected position of write iterator"); + //Update row indications + for(index_type row = smallerRow +1; row <= largerRow; ++row) { + rowIndications[row] += rowSizeDifference; + } + } + } + template ValueType SparseMatrix::getConstrainedRowSum(index_type row, storm::storage::BitVector const& constraint) const { ValueType result = storm::utility::zero(); diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 95d209b07..8e386f446 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -655,11 +655,19 @@ namespace storm { * @return True if the rows have identical entries. */ bool compareRows(index_type i1, index_type i2) const; + /*! * Finds duplicate rows in a rowgroup. */ BitVector duplicateRowsInRowgroups() const; + /** + * Swaps the two rows. + * @param row1 Index of first row + * @param row2 Index of second row + */ + void swapRows(index_type const& row1, index_type const& row2); + /*! * Selects exactly one row from each row group of this matrix and returns the resulting matrix. * diff --git a/src/utility/ModelInstantiator.h b/src/utility/ModelInstantiator.h index 939e71e2d..2b116188e 100644 --- a/src/utility/ModelInstantiator.h +++ b/src/utility/ModelInstantiator.h @@ -82,7 +82,7 @@ namespace storm { auto markovianStatesCopy = parametricModel.getMarkovianStates(); auto choiceLabelingCopy = parametricModel.getOptionalChoiceLabeling(); std::vector exitRates(parametricModel.getExitRates().size(), storm::utility::one()); - this->instantiatedModel = std::make_shared(buildDummyMatrix(parametricModel.getTransitionMatrix()), std::move(stateLabelingCopy), std::move(markovianStatesCopy), std::move(exitRates), buildDummyRewardModels(parametricModel.getRewardModels()), std::move(choiceLabelingCopy)); + this->instantiatedModel = std::make_shared(buildDummyMatrix(parametricModel.getTransitionMatrix()), std::move(stateLabelingCopy), std::move(markovianStatesCopy), std::move(exitRates), true, buildDummyRewardModels(parametricModel.getRewardModels()), std::move(choiceLabelingCopy)); initializeVectorMapping(this->instantiatedModel->getExitRates(), this->functions, this->vectorMapping, parametricModel.getExitRates()); } diff --git a/test/functional/builder/ExplicitPrismModelBuilderTest.cpp b/test/functional/builder/ExplicitPrismModelBuilderTest.cpp index ac2b1f6fb..66906c2c5 100644 --- a/test/functional/builder/ExplicitPrismModelBuilderTest.cpp +++ b/test/functional/builder/ExplicitPrismModelBuilderTest.cpp @@ -1,6 +1,7 @@ #include "gtest/gtest.h" #include "storm-config.h" #include "src/models/sparse/StandardRewardModel.h" +#include "src/models/sparse/MarkovAutomaton.h" #include "src/settings/SettingMemento.h" #include "src/parser/PrismParser.h" #include "src/builder/ExplicitModelBuilder.h" @@ -99,6 +100,31 @@ TEST(ExplicitPrismModelBuilderTest, Mdp) { EXPECT_EQ(59ul, model->getNumberOfTransitions()); } +TEST(ExplicitPrismModelBuilderTest, Ma) { + storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/simple.ma"); + + std::shared_ptr> model = storm::builder::ExplicitModelBuilder(program).build(); + EXPECT_EQ(5ul, model->getNumberOfStates()); + EXPECT_EQ(8ul, model->getNumberOfTransitions()); + ASSERT_TRUE(model->isOfType(storm::models::ModelType::MarkovAutomaton)); + EXPECT_EQ(4ul, model->as>()->getMarkovianStates().getNumberOfSetBits()); + + program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/hybrid_states.ma"); + model = storm::builder::ExplicitModelBuilder(program).build(); + EXPECT_EQ(5ul, model->getNumberOfStates()); + EXPECT_EQ(13ul, model->getNumberOfTransitions()); + ASSERT_TRUE(model->isOfType(storm::models::ModelType::MarkovAutomaton)); + EXPECT_EQ(5ul, model->as>()->getMarkovianStates().getNumberOfSetBits()); + + program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/stream2.ma"); + model = storm::builder::ExplicitModelBuilder(program).build(); + EXPECT_EQ(12ul, model->getNumberOfStates()); + EXPECT_EQ(14ul, model->getNumberOfTransitions()); + ASSERT_TRUE(model->isOfType(storm::models::ModelType::MarkovAutomaton)); + EXPECT_EQ(7ul, model->as>()->getMarkovianStates().getNumberOfSetBits()); + +} + TEST(ExplicitPrismModelBuilderTest, FailComposition) { storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/system_composition.nm"); diff --git a/test/functional/builder/hybrid_states.ma b/test/functional/builder/hybrid_states.ma new file mode 100644 index 000000000..3c19d220a --- /dev/null +++ b/test/functional/builder/hybrid_states.ma @@ -0,0 +1,19 @@ + +ma + +module hybrid_states + + s : [0..4]; + + [] (s=0) -> 0.8 : (s'=0) + 0.2 : (s'=2); + [] (s=0) -> 1 : (s' = 1); + <> (s=0) -> 3 : (s' = 1); + <> (s=1) -> 9 : (s'=0) + 1 : (s'=3); + <> (s=2) -> 12 : (s'=4); + [] (s=3) -> 1 : true; + [] (s=3) -> 1 : (s'=4); + <> (s=3) -> 2 : (s'=3) + 2 : (s'=4); + [] (s=4) -> 1 : true; + <> (s=4) -> 3 : (s'=3); + +endmodule diff --git a/test/functional/builder/simple.ma b/test/functional/builder/simple.ma new file mode 100644 index 000000000..92a003555 --- /dev/null +++ b/test/functional/builder/simple.ma @@ -0,0 +1,15 @@ + +ma + +module simple + + s : [0..4]; + + + [alpha] (s=0) -> 1 : (s' = 1); + [beta] (s=0) -> 0.8 : (s'=0) + 0.2 : (s'=2); + <> (s=1) -> 9 : (s'=0) + 1 : (s'=3); + <> (s=2) -> 12 : (s'=4); + <> (s>2) -> 1 : true; + +endmodule \ No newline at end of file diff --git a/test/functional/builder/stream2.ma b/test/functional/builder/stream2.ma new file mode 100644 index 000000000..56378a490 --- /dev/null +++ b/test/functional/builder/stream2.ma @@ -0,0 +1,50 @@ + +ma + +const int N = 2; // num packages + +const double inRate = 4; +const double processingRate = 5; + +module streamingclient + + s : [0..3]; // current state: + // 0: decide whether to start + // 1: buffering + // 2: running + // 3: done + + n : [0..N]; // number of received packages + k : [0..N]; // number of processed packages + + [buffer] s=0 & n 1 : (s'=1); + [start] s=0 & k 1 : (s'=2) & (k'=k+1); + + <> s=1 -> inRate : (n'=n+1) & (s'=0); + + <> s=2 & n inRate : (n'=n+1) + processingRate : (k'=k+1); + <> s=2 & n inRate : (n'=n+1) + processingRate : (s'=0); + <> s=2 & n=N & k processingRate : (k'=k+1); + <> s=2 & n=N & k=N -> processingRate : (s'=3); + + <> s=3 -> 1 : true; +endmodule + +// All packages received and buffer empty +label "done" = (s=3); + +rewards "buffering" + s=1 : 1; +endrewards + +rewards "initialbuffering" + s=1 & k = 0: 1; +endrewards + +rewards "intermediatebuffering" + s=1 & k > 0: 1; +endrewards + +rewards "numRestarts" + [start] k > 0 : 1; +endrewards