From 1e5398c8b7bb85a85d88a27ce204d7fc2bc01d02 Mon Sep 17 00:00:00 2001 From: dehnert Date: Sat, 20 Jun 2015 20:07:45 +0200 Subject: [PATCH] LRA finally working for ctmcs Former-commit-id: 699e4714a4cc493251ad38085fc79a93bedba75f --- .../csl/SparseCtmcCslModelChecker.cpp | 194 ++++++++++++- .../csl/SparseCtmcCslModelChecker.h | 6 +- .../prctl/SparseDtmcPrctlModelChecker.cpp | 271 +----------------- .../prctl/SparseDtmcPrctlModelChecker.h | 2 - .../modules/GmmxxEquationSolverSettings.cpp | 36 +-- .../modules/GmmxxEquationSolverSettings.h | 34 +-- .../modules/NativeEquationSolverSettings.cpp | 19 +- .../modules/NativeEquationSolverSettings.h | 18 +- src/solver/GmmxxLinearEquationSolver.cpp | 24 +- src/solver/NativeLinearEquationSolver.cpp | 131 +++++---- src/solver/NativeLinearEquationSolver.h | 5 +- src/storage/SparseMatrix.cpp | 261 ++++++++++------- src/storage/SparseMatrix.h | 9 + src/utility/solver.cpp | 21 +- src/utility/solver.h | 9 + .../NativeDtmcPrctlModelCheckerTest.cpp | 9 +- 16 files changed, 548 insertions(+), 501 deletions(-) diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp index 25328e0a5..dce356de7 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp @@ -12,6 +12,8 @@ #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" +#include "src/storage/StronglyConnectedComponentDecomposition.h" + #include "src/exceptions/InvalidStateException.h" #include "src/exceptions/InvalidPropertyException.h" #include "src/exceptions/NotImplementedException.h" @@ -476,25 +478,193 @@ namespace storm { std::unique_ptr subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - // Since we use the uniformized matrix to compute the steady state probabilities, we need to build it first. - ValueType uniformizationRate = 0; - for (auto const& rate : this->getModel().getExitRateVector()) { - uniformizationRate = std::max(uniformizationRate, rate); + storm::storage::SparseMatrix probabilityMatrix = this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(computeLongRunAverageHelper(this->getModel(), probabilityMatrix, subResult.getTruthValuesVector(), &this->getModel().getExitRateVector(), qualitative, *linearEquationSolverFactory))); + } + + template + std::vector SparseCtmcCslModelChecker::computeLongRunAverageHelper(storm::models::sparse::DeterministicModel const& model, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + // If there are no goal states, we avoid the computation and directly return zero. + uint_fast64_t numOfStates = transitionMatrix.getRowCount(); + if (psiStates.empty()) { + return std::vector(numOfStates, storm::utility::zero()); } - uniformizationRate *= 1.02; - STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive."); -// storm::storage::SparseMatrix uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), uniformizationRate, this->getModel().getExitRateVector()); - -// storm::storage::SparseMatrix uniformizedMatrix = this->computeGeneratorMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); - storm::storage::SparseMatrix uniformizedMatrix = this->computeProbabilityMatrix(this->getModel().getTransitionMatrix(), this->getModel().getExitRateVector()); + // Likewise, if all bits are set, we can avoid the computation. + if (psiStates.full()) { + return std::vector(numOfStates, storm::utility::one()); + } + + // Start by decomposing the DTMC into its BSCCs. + storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(model, false, true); + + // Get some data members for convenience. + ValueType one = storm::utility::one(); + ValueType zero = storm::utility::zero(); - std::vector result = SparseDtmcPrctlModelChecker::computeLongRunAverageHelper(this->getModel(), uniformizedMatrix, subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory); + // First we check which states are in BSCCs. + storm::storage::BitVector statesInBsccs(numOfStates); + storm::storage::BitVector firstStatesInBsccs(numOfStates); + std::vector stateToBsccIndexMap(transitionMatrix.getColumnCount()); - return std::unique_ptr(new ExplicitQuantitativeCheckResult()); + for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; + + // Gather information for later use. + bool first = true; + for (auto const& state : bscc) { + statesInBsccs.set(state); + if (first) { + firstStatesInBsccs.set(state); + } + stateToBsccIndexMap[state] = currentBsccIndex; + first = false; + } + } + + firstStatesInBsccs = firstStatesInBsccs % statesInBsccs; + storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; + + // Then we construct an equation system that yields the steady state probabilities for all states in BSCCs. + storm::storage::SparseMatrix bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true); + + // Since in the fix point equation, we need to multiply the vector from the left, we convert this to a + // multiplication from the right by transposing the system. + bsccEquationSystem = bsccEquationSystem.transpose(false, true); + + // Create an auxiliary structure that makes it easy to look up the indices within the set of BSCC states. + uint_fast64_t lastIndex = 0; + uint_fast64_t currentNumberOfSetBits = 0; + std::vector indexInStatesInBsccs; + indexInStatesInBsccs.reserve(transitionMatrix.getRowCount()); + for (auto index : statesInBsccs) { + while (lastIndex <= index) { + indexInStatesInBsccs.push_back(currentNumberOfSetBits); + ++lastIndex; + } + ++currentNumberOfSetBits; + } + + // Now build the final equation system matrix. + storm::storage::SparseMatrixBuilder builder; + uint_fast64_t currentBsccIndex = 0; + for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) { + // If the current row is the first one belonging to a BSCC, we substitute it by the constraint that the + // values for states of this BSCC must sum to one. + if (firstStatesInBsccs.get(row)) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; + + for (auto const& state : bscc) { + builder.addNextValue(row, indexInStatesInBsccs[state], one); + } + ++currentBsccIndex; + } else { + // Otherwise, we copy the row, and subtract 1 from the diagonal. + for (auto& entry : bsccEquationSystem.getRow(row)) { + if (entry.getColumn() == row) { + builder.addNextValue(row, entry.getColumn(), entry.getValue() - one); + } else { + builder.addNextValue(row, entry.getColumn(), entry.getValue()); + } + } + } + } + + bsccEquationSystem = builder.build(); + + std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); + storm::utility::vector::setVectorValues(bsccEquationSystemRightSide, firstStatesInBsccs, one); + + // As an initial guess, we take the uniform distribution over all states of the containing BSCC. + std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), zero); + for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { + storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; + + for (auto const& state : bscc) { + bsccEquationSystemSolution[indexInStatesInBsccs[state]] = one / bscc.size(); + } + } + + { + std::unique_ptr> solver = linearEquationSolverFactory.create(bsccEquationSystem); + solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide); + } + + // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour. + if (exitRateVector != nullptr) { + std::vector bsccTotalValue(bsccDecomposition.size(), zero); + std::size_t i = 0; + for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) { + bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i] * (one / (*exitRateVector)[*stateIter]); + } + + i = 0; + for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) { + bsccEquationSystemSolution[i] = (bsccEquationSystemSolution[i] * (one / (*exitRateVector)[*stateIter])) / bsccTotalValue[stateToBsccIndexMap[*stateIter]]; + } + } + + // Calculate LRA Value for each BSCC from steady state distribution in BSCCs. + std::vector bsccLra(bsccDecomposition.size(), zero); + size_t i = 0; + for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) { + if (psiStates.get(*stateIter)) { + bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i]; + } + } + + std::vector rewardSolution; + if (!statesNotInBsccs.empty()) { + // Calculate LRA for states not in bsccs as expected reachability rewards. + // Target states are states in bsccs, transition reward is the lra of the bscc for each transition into a + // bscc and 0 otherwise. This corresponds to the sum of LRAs in BSCC weighted by the reachability probability + // of the BSCC. + + std::vector rewardRightSide; + rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits()); + + for (auto state : statesNotInBsccs) { + ValueType reward = zero; + for (auto entry : transitionMatrix.getRow(state)) { + if (statesInBsccs.get(entry.getColumn())) { + reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]]; + } + } + rewardRightSide.push_back(reward); + } + + storm::storage::SparseMatrix rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true); + rewardEquationSystemMatrix.convertToEquationSystem(); + + rewardSolution = std::vector(rewardEquationSystemMatrix.getColumnCount(), one); + + { + std::unique_ptr> solver = linearEquationSolverFactory.create(rewardEquationSystemMatrix); + solver->solveEquationSystem(rewardSolution, rewardRightSide); + } + } + + // Fill the result vector. + std::vector result(numOfStates); + auto rewardSolutionIter = rewardSolution.begin(); + for (size_t state = 0; state < numOfStates; ++state) { + if (statesInBsccs.get(state)) { + // Assign the value of the bscc the state is in. + result[state] = bsccLra[stateToBsccIndexMap[state]]; + } else { + STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution."); + // Take the value from the reward computation. Since the n-th state not in any bscc is the n-th + // entry in rewardSolution we can just take the next value from the iterator. + result[state] = *rewardSolutionIter; + ++rewardSolutionIter; + } + } + + return result; } + // Explicitly instantiate the model checker. template class SparseCtmcCslModelChecker; diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.h b/src/modelchecker/csl/SparseCtmcCslModelChecker.h index 35edab8cf..42a1c9a0f 100644 --- a/src/modelchecker/csl/SparseCtmcCslModelChecker.h +++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.h @@ -13,10 +13,14 @@ namespace storm { template class HybridCtmcCslModelChecker; + template + class SparseDtmcPrctlModelChecker; + template class SparseCtmcCslModelChecker : public SparsePropositionalModelChecker { public: friend class HybridCtmcCslModelChecker; + friend class SparseDtmcPrctlModelChecker; explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc const& model); explicit SparseCtmcCslModelChecker(storm::models::sparse::Ctmc const& model, std::unique_ptr>&& linearEquationSolverFactory); @@ -40,7 +44,7 @@ namespace storm { static std::vector computeUntilProbabilitiesHelper(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); std::vector computeInstantaneousRewardsHelper(double timeBound) const; std::vector computeCumulativeRewardsHelper(double timeBound) const; - + static std::vector computeLongRunAverageHelper(storm::models::sparse::DeterministicModel const& model, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, std::vector const* exitRateVector, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); /*! * Computes the matrix representing the transitions of the uniformized CTMC. * diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 7237fe765..538aa78aa 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -8,11 +8,11 @@ #include "src/utility/graph.h" #include "src/utility/solver.h" +#include "src/modelchecker/csl/SparseCtmcCslModelChecker.h" #include "src/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "src/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "src/exceptions/InvalidStateException.h" -#include "src/storage/StronglyConnectedComponentDecomposition.h" #include "src/exceptions/InvalidPropertyException.h" @@ -302,280 +302,13 @@ namespace storm { ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeReachabilityRewardsHelper(this->getModel().getTransitionMatrix(), this->getModel().getOptionalStateRewardVector(), this->getModel().getOptionalTransitionRewardMatrix(), this->getModel().getBackwardTransitions(), subResult.getTruthValuesVector(), *this->linearEquationSolverFactory, qualitative))); } - - template - std::vector SparseDtmcPrctlModelChecker::computeLongRunAverageHelper(storm::models::sparse::DeterministicModel const& model, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - // If there are no goal states, we avoid the computation and directly return zero. - uint_fast64_t numOfStates = transitionMatrix.getRowCount(); - if (psiStates.empty()) { - return std::vector(numOfStates, storm::utility::zero()); - } - - // Likewise, if all bits are set, we can avoid the computation. - if (psiStates.full()) { - return std::vector(numOfStates, storm::utility::one()); - } - - // Start by decomposing the DTMC into its BSCCs. - storm::storage::StronglyConnectedComponentDecomposition bsccDecomposition(model, false, true); - - // Get some data members for convenience. - ValueType one = storm::utility::one(); - ValueType zero = storm::utility::zero(); - - // First we check which states are in BSCCs. We use this later to speed up reachability analysis - storm::storage::BitVector statesInBsccs(numOfStates); - storm::storage::BitVector statesInBsccsWithoutFirst(numOfStates); - - std::vector stateToBsccIndexMap(transitionMatrix.getColumnCount()); - - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - - // Gather information for later use. - bool first = true; - for (auto const& state : bscc) { - statesInBsccs.set(state); - if (!first) { - statesInBsccsWithoutFirst.set(state); - } - stateToBsccIndexMap[state] = currentBsccIndex; - first = false; - } - } - - storm::storage::BitVector statesNotInBsccs = ~statesInBsccs; - - std::cout << transitionMatrix << std::endl; - - // Calculate steady state distribution for all BSCCs by calculating an eigenvector for the eigenvalue 1 of - // the transposed transition matrix for the BSCCs. - storm::storage::SparseMatrix bsccEquationSystem = transitionMatrix.getSubmatrix(false, statesInBsccs, statesInBsccs, true); - - std::cout << bsccEquationSystem << std::endl; - - bsccEquationSystem = bsccEquationSystem.transpose(false, true); - - std::cout << bsccEquationSystem << std::endl; - - // Subtract identity matrix. - for (uint_fast64_t row = 0; row < bsccEquationSystem.getRowCount(); ++row) { - for (auto& entry : bsccEquationSystem.getRow(row)) { - if (entry.getColumn() == row) { - entry.setValue(entry.getValue() - one); - } - } - } - - std::cout << bsccEquationSystem << std::endl; - - std::cout << statesInBsccsWithoutFirst << " // " << statesInBsccs << std::endl; - bsccEquationSystem = bsccEquationSystem.getSubmatrix(false, statesInBsccsWithoutFirst, statesInBsccs, false); - - std::cout << bsccEquationSystem << std::endl; - - // For each BSCC, add a row to the matrix that expresses that the sum over all entries in this BSCC needs to be one. - storm::storage::SparseMatrixBuilder builder(std::move(bsccEquationSystem)); - typename storm::storage::SparseMatrixBuilder::index_type row = builder.getLastRow() + 1; - - for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - - for (auto const& state : bscc) { - builder.addNextValue(row, state, one); - } - ++row; - } - bsccEquationSystem = builder.build(); - - std::cout << bsccEquationSystem << std::endl; - - std::vector bsccEquationSystemRightSide(bsccEquationSystem.getColumnCount(), zero); - bsccEquationSystemRightSide.back() = one; - std::vector bsccEquationSystemSolution(bsccEquationSystem.getColumnCount(), one); - { - std::unique_ptr> solver = linearEquationSolverFactory.create(bsccEquationSystem); - solver->solveEquationSystem(bsccEquationSystemSolution, bsccEquationSystemRightSide); - } - - ValueType sum = zero; - for (auto const& elem : bsccEquationSystemSolution) { - std::cout << "sol " << elem << std::endl; - sum += elem; - } - std::cout << "sum: " << sum << std::endl; - std::cout << "in " << bsccDecomposition.size() << "bsccs" << std::endl; - - // Calculate LRA Value for each BSCC from steady state distribution in BSCCs. - // We have to scale the results, as the probabilities for each BSCC have to sum up to one, which they don't - // necessarily do in the solution of the equation system. - std::vector bsccLra(bsccDecomposition.size(), zero); - std::vector bsccTotalValue(bsccDecomposition.size(), zero); - size_t i = 0; - for (auto stateIter = statesInBsccs.begin(); stateIter != statesInBsccs.end(); ++i, ++stateIter) { - if (psiStates.get(*stateIter)) { - bsccLra[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i]; - } - bsccTotalValue[stateToBsccIndexMap[*stateIter]] += bsccEquationSystemSolution[i]; - } - for (i = 0; i < bsccLra.size(); ++i) { - bsccLra[i] /= bsccTotalValue[i]; - } - - std::vector rewardSolution; - if (!statesNotInBsccs.empty()) { - // Calculate LRA for states not in bsccs as expected reachability rewards. - // Target states are states in bsccs, transition reward is the lra of the bscc for each transition into a - // bscc and 0 otherwise. This corresponds to the sum of LRAs in BSCC weighted by the reachability probability - // of the BSCC. - - std::vector rewardRightSide; - rewardRightSide.reserve(statesNotInBsccs.getNumberOfSetBits()); - - for (auto state : statesNotInBsccs) { - ValueType reward = zero; - for (auto entry : transitionMatrix.getRow(state)) { - if (statesInBsccs.get(entry.getColumn())) { - reward += entry.getValue() * bsccLra[stateToBsccIndexMap[entry.getColumn()]]; - } - } - rewardRightSide.push_back(reward); - } - - storm::storage::SparseMatrix rewardEquationSystemMatrix = transitionMatrix.getSubmatrix(false, statesNotInBsccs, statesNotInBsccs, true); - rewardEquationSystemMatrix.convertToEquationSystem(); - - rewardSolution = std::vector(rewardEquationSystemMatrix.getColumnCount(), one); - - { - std::unique_ptr> solver = linearEquationSolverFactory.create(rewardEquationSystemMatrix); - solver->solveEquationSystem(rewardSolution, rewardRightSide); - } - } - - // Fill the result vector. - std::vector result(numOfStates); - auto rewardSolutionIter = rewardSolution.begin(); - for (size_t state = 0; state < numOfStates; ++state) { - if (statesInBsccs.get(state)) { - // Assign the value of the bscc the state is in. - result[state] = bsccLra[stateToBsccIndexMap[state]]; - } else { - STORM_LOG_ASSERT(rewardSolutionIter != rewardSolution.end(), "Too few elements in solution."); - // Take the value from the reward computation. Since the n-th state not in any bscc is the n-th - // entry in rewardSolution we can just take the next value from the iterator. - result[state] = *rewardSolutionIter; - ++rewardSolutionIter; - } - } - - return result; - - //old implementeation - - //now we calculate the long run average for each BSCC in isolation and compute the weighted contribution of the BSCC to the LRA value of all states - //for (uint_fast64_t currentBsccIndex = 0; currentBsccIndex < bsccDecomposition.size(); ++currentBsccIndex) { - // storm::storage::StronglyConnectedComponent const& bscc = bsccDecomposition[currentBsccIndex]; - - // storm::storage::BitVector statesInThisBscc(numOfStates); - // for (auto const& state : bscc) { - // statesInThisBscc.set(state); - // } - - // //ValueType lraForThisBscc = this->computeLraForBSCC(transitionMatrix, psiStates, bscc); - // ValueType lraForThisBscc = bsccLra[currentBsccIndex]; - - // //the LRA value of a BSCC contributes to the LRA value of a state with the probability of reaching that BSCC from that state - // //thus we add Prob(Eventually statesInThisBscc) * lraForThisBscc to the result vector - - // //the reachability probabilities will be zero in other BSCCs, thus we can set the left operand of the until formula to statesNotInBsccs as an optimization - // std::vector reachProbThisBscc = this->computeUntilProbabilitiesHelper(statesNotInBsccs, statesInThisBscc, false); - // - // storm::utility::vector::scaleVectorInPlace(reachProbThisBscc, lraForThisBscc); - // storm::utility::vector::addVectorsInPlace(result, reachProbThisBscc); - //} - - //return result; - } template std::unique_ptr SparseDtmcPrctlModelChecker::computeLongRunAverage(storm::logic::StateFormula const& stateFormula, bool qualitative, boost::optional const& optimalityType) { std::unique_ptr subResultPointer = this->check(stateFormula); ExplicitQualitativeCheckResult const& subResult = subResultPointer->asExplicitQualitativeCheckResult(); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(this->computeLongRunAverageHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), qualitative, *linearEquationSolverFactory))); - } - - - template - ValueType SparseDtmcPrctlModelChecker::computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { - //if size is one we can immediately derive the result - if (bscc.size() == 1){ - if (psiStates.get(*(bscc.begin()))) { - return storm::utility::one(); - } else{ - return storm::utility::zero(); - } - } - - storm::storage::BitVector subsystem = storm::storage::BitVector(transitionMatrix.getRowCount()); - subsystem.set(bscc.begin(), bscc.end()); - - //we now have to solve ((P')^t - I) * x = 0, where P' is the submatrix of transitionMatrix, - // ^t means transose, and I is the identity matrix. - - storm::storage::SparseMatrix subsystemMatrix = transitionMatrix.getSubmatrix(false, subsystem, subsystem, true); - subsystemMatrix = subsystemMatrix.transpose(); - - // subtract 1 on the diagonal and at the same time add a row with all ones to enforce that the result of the equation system is unique - storm::storage::SparseMatrixBuilder equationSystemBuilder(subsystemMatrix.getRowCount() + 1, subsystemMatrix.getColumnCount(), subsystemMatrix.getEntryCount() + subsystemMatrix.getColumnCount()); - ValueType one = storm::utility::one(); - ValueType zero = storm::utility::zero(); - bool foundDiagonalElement = false; - for (uint_fast64_t row = 0; row < subsystemMatrix.getRowCount(); ++row) { - for (auto& entry : subsystemMatrix.getRowGroup(row)) { - if (entry.getColumn() == row) { - equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue() - one); - foundDiagonalElement = true; - } else { - equationSystemBuilder.addNextValue(row, entry.getColumn(), entry.getValue()); - } - } - - // Throw an exception if a row did not have an element on the diagonal. - STORM_LOG_THROW(foundDiagonalElement, storm::exceptions::InvalidOperationException, "Internal Error, no diagonal entry found."); - } - for (uint_fast64_t column = 0; column + 1 < subsystemMatrix.getColumnCount(); ++column) { - equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), column, one); - } - equationSystemBuilder.addNextValue(subsystemMatrix.getRowCount(), subsystemMatrix.getColumnCount() - 1, zero); - subsystemMatrix = equationSystemBuilder.build(); - - // create x and b for the equation system. setting the last entry of b to one enforces the sum over the unique solution vector is one - std::vector steadyStateDist(subsystemMatrix.getRowCount(), zero); - std::vector b(subsystemMatrix.getRowCount(), zero); - b[subsystemMatrix.getRowCount() - 1] = one; - - - { - auto solver = linearEquationSolverFactory.create(subsystemMatrix); - solver->solveEquationSystem(steadyStateDist, b); - } - - //remove the last entry of the vector as it was just there to enforce the unique solution - steadyStateDist.pop_back(); - - //calculate the fraction we spend in psi states on the long run - std::vector steadyStateForPsiStates(transitionMatrix.getRowCount() - 1, zero); - storm::utility::vector::setVectorValues(steadyStateForPsiStates, psiStates & subsystem, steadyStateDist); - - ValueType result = zero; - - for (auto value : steadyStateForPsiStates) { - result += value; - } - - return result; + return std::unique_ptr(new ExplicitQuantitativeCheckResult(SparseCtmcCslModelChecker::computeLongRunAverageHelper(this->getModel(), this->getModel().getTransitionMatrix(), subResult.getTruthValuesVector(), nullptr, qualitative, *linearEquationSolverFactory))); } template diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h index 212b99588..5409091fd 100644 --- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h +++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h @@ -43,8 +43,6 @@ namespace storm { static std::vector computeReachabilityRewardsHelper(storm::storage::SparseMatrix const& transitionMatrix, boost::optional> const& stateRewardVector, boost::optional> const& transitionRewardMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& targetStates, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, bool qualitative); static std::vector computeLongRunAverageHelper(storm::models::sparse::DeterministicModel const& model, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& psiStates, bool qualitative, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - static ValueType computeLraForBSCC(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& goalStates, storm::storage::StronglyConnectedComponent const& bscc, storm::utility::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); - // An object that is used for retrieving linear equation solvers. std::unique_ptr> linearEquationSolverFactory; }; diff --git a/src/settings/modules/GmmxxEquationSolverSettings.cpp b/src/settings/modules/GmmxxEquationSolverSettings.cpp index 4b9f27f9e..359836968 100644 --- a/src/settings/modules/GmmxxEquationSolverSettings.cpp +++ b/src/settings/modules/GmmxxEquationSolverSettings.cpp @@ -7,7 +7,7 @@ namespace storm { namespace modules { const std::string GmmxxEquationSolverSettings::moduleName = "gmm++"; - const std::string GmmxxEquationSolverSettings::techniqueOptionName = "tech"; + const std::string GmmxxEquationSolverSettings::techniqueOptionName = "method"; const std::string GmmxxEquationSolverSettings::preconditionOptionName = "precond"; const std::string GmmxxEquationSolverSettings::restartOptionName = "restart"; const std::string GmmxxEquationSolverSettings::maximalIterationsOptionName = "maxiter"; @@ -32,38 +32,38 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build()); } - bool GmmxxEquationSolverSettings::isLinearEquationSystemTechniqueSet() const { + bool GmmxxEquationSolverSettings::isLinearEquationSystemMethodSet() const { return this->getOption(techniqueOptionName).getHasOptionBeenSet(); } - GmmxxEquationSolverSettings::LinearEquationTechnique GmmxxEquationSolverSettings::getLinearEquationSystemTechnique() const { + GmmxxEquationSolverSettings::LinearEquationMethod GmmxxEquationSolverSettings::getLinearEquationSystemMethod() const { std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString(); if (linearEquationSystemTechniqueAsString == "bicgstab") { - return GmmxxEquationSolverSettings::LinearEquationTechnique::Bicgstab; + return GmmxxEquationSolverSettings::LinearEquationMethod::Bicgstab; } else if (linearEquationSystemTechniqueAsString == "qmr") { - return GmmxxEquationSolverSettings::LinearEquationTechnique::Qmr; + return GmmxxEquationSolverSettings::LinearEquationMethod::Qmr; } else if (linearEquationSystemTechniqueAsString == "gmres") { - return GmmxxEquationSolverSettings::LinearEquationTechnique::Gmres; + return GmmxxEquationSolverSettings::LinearEquationMethod::Gmres; } else if (linearEquationSystemTechniqueAsString == "jacobi") { - return GmmxxEquationSolverSettings::LinearEquationTechnique::Jacobi; + return GmmxxEquationSolverSettings::LinearEquationMethod::Jacobi; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected."); } - bool GmmxxEquationSolverSettings::isPreconditioningTechniqueSet() const { + bool GmmxxEquationSolverSettings::isPreconditioningMethodSet() const { return this->getOption(preconditionOptionName).getHasOptionBeenSet(); } - GmmxxEquationSolverSettings::PreconditioningTechnique GmmxxEquationSolverSettings::getPreconditioningTechnique() const { - std::string preconditioningTechniqueAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString(); - if (preconditioningTechniqueAsString == "ilu") { - return GmmxxEquationSolverSettings::PreconditioningTechnique::Ilu; - } else if (preconditioningTechniqueAsString == "diagonal") { - return GmmxxEquationSolverSettings::PreconditioningTechnique::Diagonal; - } else if (preconditioningTechniqueAsString == "none") { - return GmmxxEquationSolverSettings::PreconditioningTechnique::None; + GmmxxEquationSolverSettings::PreconditioningMethod GmmxxEquationSolverSettings::getPreconditioningMethod() const { + std::string PreconditioningMethodAsString = this->getOption(preconditionOptionName).getArgumentByName("name").getValueAsString(); + if (PreconditioningMethodAsString == "ilu") { + return GmmxxEquationSolverSettings::PreconditioningMethod::Ilu; + } else if (PreconditioningMethodAsString == "diagonal") { + return GmmxxEquationSolverSettings::PreconditioningMethod::Diagonal; + } else if (PreconditioningMethodAsString == "none") { + return GmmxxEquationSolverSettings::PreconditioningMethod::None; } - STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << preconditioningTechniqueAsString << "' selected."); + STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown preconditioning technique '" << PreconditioningMethodAsString << "' selected."); } bool GmmxxEquationSolverSettings::isRestartIterationCountSet() const { @@ -99,7 +99,7 @@ namespace storm { } bool GmmxxEquationSolverSettings::check() const { - bool optionsSet = isLinearEquationSystemTechniqueSet() || isPreconditioningTechniqueSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet(); + bool optionsSet = isLinearEquationSystemMethodSet() || isPreconditioningMethodSet() || isRestartIterationCountSet() | isMaximalIterationCountSet() || isPrecisionSet() || isConvergenceCriterionSet(); STORM_LOG_WARN_COND(storm::settings::generalSettings().getEquationSolver() == storm::settings::modules::GeneralSettings::EquationSolver::Gmmxx || !optionsSet, "gmm++ is not selected as the equation solver, so setting options for gmm++ has no effect."); diff --git a/src/settings/modules/GmmxxEquationSolverSettings.h b/src/settings/modules/GmmxxEquationSolverSettings.h index 07bc90174..8cb0761a8 100644 --- a/src/settings/modules/GmmxxEquationSolverSettings.h +++ b/src/settings/modules/GmmxxEquationSolverSettings.h @@ -12,11 +12,11 @@ namespace storm { */ class GmmxxEquationSolverSettings : public ModuleSettings { public: - // An enumeration of all available techniques for solving linear equations. - enum class LinearEquationTechnique { Bicgstab, Qmr, Gmres, Jacobi }; + // An enumeration of all available methods for solving linear equations. + enum class LinearEquationMethod { Bicgstab, Qmr, Gmres, Jacobi }; - // An enumeration of all available preconditioning techniques. - enum class PreconditioningTechnique { Ilu, Diagonal, None }; + // An enumeration of all available preconditioning methods. + enum class PreconditioningMethod { Ilu, Diagonal, None }; // An enumeration of all available convergence criteria. enum class ConvergenceCriterion { Absolute, Relative }; @@ -29,32 +29,32 @@ namespace storm { GmmxxEquationSolverSettings(storm::settings::SettingsManager& settingsManager); /*! - * Retrieves whether the linear equation system technique has been set. + * Retrieves whether the linear equation system method has been set. * - * @return True iff the linear equation system technique has been set. + * @return True iff the linear equation system method has been set. */ - bool isLinearEquationSystemTechniqueSet() const; + bool isLinearEquationSystemMethodSet() const; /*! - * Retrieves the technique that is to be used for solving systems of linear equations. + * Retrieves the method that is to be used for solving systems of linear equations. * - * @return The technique to use. + * @return The method to use. */ - LinearEquationTechnique getLinearEquationSystemTechnique() const; + LinearEquationMethod getLinearEquationSystemMethod() const; /*! - * Retrieves whether the preconditioning technique has been set. + * Retrieves whether the preconditioning method has been set. * - * @return True iff the preconditioning technique has been set. + * @return True iff the preconditioning method has been set. */ - bool isPreconditioningTechniqueSet() const; + bool isPreconditioningMethodSet() const; /*! - * Retrieves the technique that is to be used for preconditioning solving systems of linear equations. + * Retrieves the method that is to be used for preconditioning solving systems of linear equations. * - * @return The technique to use. + * @return The method to use. */ - PreconditioningTechnique getPreconditioningTechnique() const; + PreconditioningMethod getPreconditioningMethod() const; /*! * Retrieves whether the restart iteration count has been set. @@ -64,7 +64,7 @@ namespace storm { bool isRestartIterationCountSet() const; /*! - * Retrieves the number of iterations after which restarted techniques are to be restarted. + * Retrieves the number of iterations after which restarted methods are to be restarted. * * @return The number of iterations after which to restart. */ diff --git a/src/settings/modules/NativeEquationSolverSettings.cpp b/src/settings/modules/NativeEquationSolverSettings.cpp index 76af99b3e..482713ed8 100644 --- a/src/settings/modules/NativeEquationSolverSettings.cpp +++ b/src/settings/modules/NativeEquationSolverSettings.cpp @@ -7,20 +7,23 @@ namespace storm { namespace modules { const std::string NativeEquationSolverSettings::moduleName = "native"; - const std::string NativeEquationSolverSettings::techniqueOptionName = "tech"; + const std::string NativeEquationSolverSettings::techniqueOptionName = "method"; + const std::string NativeEquationSolverSettings::omegaOptionName = "soromega"; const std::string NativeEquationSolverSettings::maximalIterationsOptionName = "maxiter"; const std::string NativeEquationSolverSettings::maximalIterationsOptionShortName = "i"; const std::string NativeEquationSolverSettings::precisionOptionName = "precision"; const std::string NativeEquationSolverSettings::absoluteOptionName = "absolute"; NativeEquationSolverSettings::NativeEquationSolverSettings(storm::settings::SettingsManager& settingsManager) : ModuleSettings(settingsManager, moduleName) { - std::vector methods = { "jacobi" }; + std::vector methods = { "jacobi", "gaussseidel", "sor" }; this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the native engine. Available are: { jacobi }.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidationFunctionString(storm::settings::ArgumentValidators::stringInListValidator(methods)).setDefaultValueString("jacobi").build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, omegaOptionName, false, "The omega used for SOR.").addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The value of the SOR parameter.").setDefaultValueDouble(0.9).addValidationFunctionDouble(storm::settings::ArgumentValidators::doubleRangeValidatorExcluding(0.0, 1.0)).build()).build()); + this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").build()); } @@ -28,10 +31,14 @@ namespace storm { return this->getOption(techniqueOptionName).getHasOptionBeenSet(); } - NativeEquationSolverSettings::LinearEquationTechnique NativeEquationSolverSettings::getLinearEquationSystemTechnique() const { + NativeEquationSolverSettings::LinearEquationMethod NativeEquationSolverSettings::getLinearEquationSystemMethod() const { std::string linearEquationSystemTechniqueAsString = this->getOption(techniqueOptionName).getArgumentByName("name").getValueAsString(); if (linearEquationSystemTechniqueAsString == "jacobi") { - return NativeEquationSolverSettings::LinearEquationTechnique::Jacobi; + return NativeEquationSolverSettings::LinearEquationMethod::Jacobi; + } else if (linearEquationSystemTechniqueAsString == "gaussseidel") { + return NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel; + } else if (linearEquationSystemTechniqueAsString == "sor") { + return NativeEquationSolverSettings::LinearEquationMethod::SOR; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected."); } @@ -51,6 +58,10 @@ namespace storm { double NativeEquationSolverSettings::getPrecision() const { return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble(); } + + double NativeEquationSolverSettings::getOmega() const { + return this->getOption(omegaOptionName).getArgumentByName("value").getValueAsDouble(); + } bool NativeEquationSolverSettings::isConvergenceCriterionSet() const { return this->getOption(absoluteOptionName).getHasOptionBeenSet(); diff --git a/src/settings/modules/NativeEquationSolverSettings.h b/src/settings/modules/NativeEquationSolverSettings.h index 9d6e43ecd..01c0fd7c9 100644 --- a/src/settings/modules/NativeEquationSolverSettings.h +++ b/src/settings/modules/NativeEquationSolverSettings.h @@ -12,8 +12,8 @@ namespace storm { */ class NativeEquationSolverSettings : public ModuleSettings { public: - // An enumeration of all available techniques for solving linear equations. - enum class LinearEquationTechnique { Jacobi }; + // An enumeration of all available methods for solving linear equations. + enum class LinearEquationMethod { Jacobi, GaussSeidel, SOR }; // An enumeration of all available convergence criteria. enum class ConvergenceCriterion { Absolute, Relative }; @@ -33,11 +33,11 @@ namespace storm { bool isLinearEquationSystemTechniqueSet() const; /*! - * Retrieves the technique that is to be used for solving systems of linear equations. + * Retrieves the method that is to be used for solving systems of linear equations. * - * @return The technique to use. + * @return The method to use. */ - LinearEquationTechnique getLinearEquationSystemTechnique() const; + LinearEquationMethod getLinearEquationSystemMethod() const; /*! * Retrieves whether the maximal iteration count has been set. @@ -67,6 +67,13 @@ namespace storm { */ double getPrecision() const; + /*! + * Retrieves the value of omega to be used for the SOR method. + * + * @return The value of omega to be used. + */ + double getOmega() const; + /*! * Retrieves whether the convergence criterion has been set. * @@ -89,6 +96,7 @@ namespace storm { private: // Define the string names of the options as constants. static const std::string techniqueOptionName; + static const std::string omegaOptionName; static const std::string maximalIterationsOptionName; static const std::string maximalIterationsOptionShortName; static const std::string precisionOptionName; diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp index cbf2eafc8..bd43c3e88 100644 --- a/src/solver/GmmxxLinearEquationSolver.cpp +++ b/src/solver/GmmxxLinearEquationSolver.cpp @@ -31,24 +31,24 @@ namespace storm { restart = settings.getRestartIterationCount(); // Determine the method to be used. - storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique methodAsSetting = settings.getLinearEquationSystemTechnique(); - if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Bicgstab) { + storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod methodAsSetting = settings.getLinearEquationSystemMethod(); + if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Bicgstab) { method = SolutionMethod::Bicgstab; - } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Qmr) { + } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Qmr) { method = SolutionMethod::Qmr; - } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Gmres) { + } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Gmres) { method = SolutionMethod::Gmres; - } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationTechnique::Jacobi) { + } else if (methodAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::LinearEquationMethod::Jacobi) { method = SolutionMethod::Jacobi; } // Check which preconditioner to use. - storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique preconditionAsSetting = settings.getPreconditioningTechnique(); - if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::Ilu) { + storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod preconditionAsSetting = settings.getPreconditioningMethod(); + if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::Ilu) { preconditioner = Preconditioner::Ilu; - } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::Diagonal) { + } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::Diagonal) { preconditioner = Preconditioner::Diagonal; - } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningTechnique::None) { + } else if (preconditionAsSetting == storm::settings::modules::GmmxxEquationSolverSettings::PreconditioningMethod::None) { preconditioner = Preconditioner::None; } } @@ -149,12 +149,6 @@ namespace storm { uint_fast64_t GmmxxLinearEquationSolver::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector* multiplyResult) const { // Get a Jacobi decomposition of the matrix A. std::pair, std::vector> jacobiDecomposition = A.getJacobiDecomposition(); - std::cout << "LU" << std::endl; - std::cout << jacobiDecomposition.first << std::endl; - for (auto const& elem : jacobiDecomposition.second) { - std::cout << elem << std::endl; - } - std::cout << "----" << std::endl; // Convert the LU matrix to gmm++'s format. std::unique_ptr> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix(std::move(jacobiDecomposition.first)); diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp index a2480ddc5..a389853c6 100644 --- a/src/solver/NativeLinearEquationSolver.cpp +++ b/src/solver/NativeLinearEquationSolver.cpp @@ -15,7 +15,7 @@ namespace storm { } template - NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A) : A(A) { + NativeLinearEquationSolver::NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method) : A(A), method(method) { // Get the settings object to customize linear solving. storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::nativeEquationSolverSettings(); @@ -23,61 +23,94 @@ namespace storm { maximalNumberOfIterations = settings.getMaximalIterationCount(); precision = settings.getPrecision(); relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative; - - // Determine the method to be used. - storm::settings::modules::NativeEquationSolverSettings::LinearEquationTechnique methodAsSetting = settings.getLinearEquationSystemTechnique(); - if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationTechnique::Jacobi) { - method = SolutionMethod::Jacobi; - } } template void NativeLinearEquationSolver::solveEquationSystem(std::vector& x, std::vector const& b, std::vector* multiplyResult) const { - // Get a Jacobi decomposition of the matrix A. - 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 precision calculation. - std::vector tmpX(x.size()); - - // Set up additional environment variables. - uint_fast64_t iterationCount = 0; - bool converged = false; - - while (!converged && iterationCount < maximalNumberOfIterations) { - // 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); + if (method == SolutionMethod::SOR || method == SolutionMethod::GaussSeidel) { + // Define the omega used for SOR. + ValueType omega = method == SolutionMethod::SOR ? storm::settings::nativeEquationSolverSettings().getOmega() : storm::utility::one(); - // Swap the two pointers as a preparation for the next iteration. - std::swap(nextX, currentX); + // 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; + } - // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, static_cast(precision), relative); + // Set up additional environment variables. + uint_fast64_t iterationCount = 0; + bool converged = false; - // Increase iteration count so we can abort if convergence is too slow. - ++iterationCount; - } - - // 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) { - std::swap(x, *currentX); - } - - // If the vector for the temporary multiplication result was not provided, we need to delete it. - if (!multiplyResultProvided) { - delete copyX; + while (!converged && iterationCount < maximalNumberOfIterations) { + A.performSuccessiveOverRelaxationStep(omega, x, b); + + // Now check if the process already converged within our precision. + converged = storm::utility::vector::equalModuloPrecision(x, *tmpX, static_cast(precision), relative); + + // If we did not yet converge, we need to copy the contents of x to *tmpX. + if (!converged) { + *tmpX = 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(); + + // 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()); + + // Set up additional environment variables. + uint_fast64_t iterationCount = 0; + bool converged = false; + + while (!converged && iterationCount < maximalNumberOfIterations) { + // 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); + + // Now check if the process already converged within our precision. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, static_cast(precision), relative); + + // Increase iteration count so we can abort if convergence is too slow. + ++iterationCount; + } + + // 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) { + std::swap(x, *currentX); + } + + // If the vector for the temporary multiplication result was not provided, we need to delete it. + if (!multiplyResultProvided) { + delete copyX; + } } } diff --git a/src/solver/NativeLinearEquationSolver.h b/src/solver/NativeLinearEquationSolver.h index f8f5a1378..aa20b1e71 100644 --- a/src/solver/NativeLinearEquationSolver.h +++ b/src/solver/NativeLinearEquationSolver.h @@ -14,15 +14,16 @@ namespace storm { public: // An enumeration specifying the available solution methods. enum class SolutionMethod { - Jacobi + Jacobi, GaussSeidel, SOR }; /*! * Constructs a linear equation solver with parameters being set according to the settings object. * * @param A The matrix defining the coefficients of the linear equation system. + * @param method The method to use for solving linear equations. */ - NativeLinearEquationSolver(storm::storage::SparseMatrix const& A); + NativeLinearEquationSolver(storm::storage::SparseMatrix const& A, SolutionMethod method = SolutionMethod::Jacobi); /*! * Constructs a linear equation solver with the given parameters. diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 2db8b576b..3dd8910dc 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -50,7 +50,7 @@ namespace storm { void MatrixEntry::setValue(ValueType const& value) { this->entry.second = value; } - + template std::pair const& MatrixEntry::getColumnValuePair() const { return this->entry; @@ -102,7 +102,7 @@ namespace storm { if (row < lastRow) { throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrixBuilder::addNextValue: adding an element in row " << row << ", but an element in row " << lastRow << " has already been added."; } - + // Check that we did not move backwards wrt. to column. if (row == lastRow && column < lastColumn) { throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrixBuilder::addNextValue: adding an element in column " << column << " in row " << row << ", but an element in column " << lastColumn << " has already been added in that row."; @@ -117,14 +117,14 @@ namespace storm { lastRow = row; } - + lastColumn = column; // Finally, set the element and increase the current size. columnsAndValues.emplace_back(column, value); highestColumn = std::max(highestColumn, column); ++currentEntryCount; - + // In case we did not expect this value, we throw an exception. if (forceInitialDimensions) { STORM_LOG_THROW(!initialRowCountSet || lastRow < initialRowCount, storm::exceptions::OutOfRangeException, "Cannot insert value at illegal row " << lastRow << "."); @@ -159,14 +159,14 @@ namespace storm { // as now the indices of row i are always between rowIndications[i] and rowIndications[i + 1], also for // the first and last row. rowIndications.push_back(currentEntryCount); - + uint_fast64_t columnCount = highestColumn + 1; if (initialColumnCountSet && forceInitialDimensions) { STORM_LOG_THROW(columnCount <= initialColumnCount, storm::exceptions::InvalidStateException, "Expected not more than " << initialColumnCount << " columns, but got " << columnCount << "."); columnCount = std::max(columnCount, initialColumnCount); } columnCount = std::max(columnCount, overriddenColumnCount); - + uint_fast64_t entryCount = currentEntryCount; if (initialEntryCountSet && forceInitialDimensions) { STORM_LOG_THROW(entryCount == initialEntryCount, storm::exceptions::InvalidStateException, "Expected " << initialEntryCount << " entries, but got " << entryCount << "."); @@ -184,12 +184,12 @@ namespace storm { rowGroupCount = std::max(rowGroupCount, initialRowGroupCount); } rowGroupCount = std::max(rowGroupCount, overriddenRowGroupCount); - + for (index_type i = currentRowGroup; i <= rowGroupCount; ++i) { rowGroupIndices.push_back(rowCount); } } - + return SparseMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), hasCustomRowGrouping); } @@ -270,14 +270,14 @@ namespace storm { template SparseMatrix::SparseMatrix(index_type columnCount, std::vector const& rowIndications, std::vector> const& columnsAndValues, std::vector const& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(rowGroupIndices) { - this->updateNonzeroEntryCount(); + this->updateNonzeroEntryCount(); } template SparseMatrix::SparseMatrix(index_type columnCount, std::vector&& rowIndications, std::vector>&& columnsAndValues, std::vector&& rowGroupIndices, bool nontrivialRowGrouping) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), nontrivialRowGrouping(nontrivialRowGrouping), rowGroupIndices(std::move(rowGroupIndices)) { - this->updateNonzeroEntryCount(); + this->updateNonzeroEntryCount(); } - + template SparseMatrix& SparseMatrix::operator=(SparseMatrix const& other) { // Only perform assignment if source and target are not the same. @@ -321,7 +321,7 @@ namespace storm { } bool equalityResult = true; - + equalityResult &= this->getRowCount() == other.getRowCount(); if (!equalityResult) { return false; @@ -379,35 +379,35 @@ namespace storm { return entryCount; } - template - uint_fast64_t SparseMatrix::getRowGroupEntryCount(uint_fast64_t const group) const { - uint_fast64_t result = 0; - for (uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group + 1]; ++row) { - result += (this->rowIndications[row + 1] - this->rowIndications[row]); - } - return result; - } - + template + uint_fast64_t SparseMatrix::getRowGroupEntryCount(uint_fast64_t const group) const { + uint_fast64_t result = 0; + for (uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group + 1]; ++row) { + result += (this->rowIndications[row + 1] - this->rowIndications[row]); + } + return result; + } + template typename SparseMatrix::index_type SparseMatrix::getNonzeroEntryCount() const { return nonzeroEntryCount; - } - - template - void SparseMatrix::updateNonzeroEntryCount() const { - //SparseMatrix* self = const_cast*>(this); - this->nonzeroEntryCount = 0; - for (auto const& element : *this) { - if (element.getValue() != storm::utility::zero()) { - ++this->nonzeroEntryCount; - } - } - } - - template - void SparseMatrix::updateNonzeroEntryCount(std::make_signed::type difference) { - this->nonzeroEntryCount += difference; - } + } + + template + void SparseMatrix::updateNonzeroEntryCount() const { + //SparseMatrix* self = const_cast*>(this); + this->nonzeroEntryCount = 0; + for (auto const& element : *this) { + if (element.getValue() != storm::utility::zero()) { + ++this->nonzeroEntryCount; + } + } + } + + template + void SparseMatrix::updateNonzeroEntryCount(std::make_signed::type difference) { + this->nonzeroEntryCount += difference; + } template typename SparseMatrix::index_type SparseMatrix::getRowGroupCount() const { @@ -423,7 +423,7 @@ namespace storm { std::vector::index_type> const& SparseMatrix::getRowGroupIndices() const { return rowGroupIndices; } - + template void SparseMatrix::makeRowsAbsorbing(storm::storage::BitVector const& rows) { for (auto row : rows) { @@ -457,7 +457,7 @@ namespace storm { columnValuePtr->setValue(storm::utility::one()); ++columnValuePtr; for (; columnValuePtr != columnValuePtrEnd; ++columnValuePtr) { - ++this->nonzeroEntryCount; + ++this->nonzeroEntryCount; columnValuePtr->setColumn(0); columnValuePtr->setValue(storm::utility::zero()); } @@ -510,33 +510,42 @@ namespace storm { return getSubmatrix(rowConstraint, columnConstraint, fakeRowGroupIndices, insertDiagonalElements); } } - + template SparseMatrix SparseMatrix::getSubmatrix(storm::storage::BitVector const& rowGroupConstraint, storm::storage::BitVector const& columnConstraint, std::vector const& rowGroupIndices, bool insertDiagonalEntries) const { uint_fast64_t submatrixColumnCount = columnConstraint.getNumberOfSetBits(); // Start by creating a temporary vector that stores for each index whose bit is set to true the number of // bits that were set before that particular index. - std::vector bitsSetBeforeIndex; - bitsSetBeforeIndex.reserve(columnCount); + std::vector columnBitsSetBeforeIndex; + columnBitsSetBeforeIndex.reserve(columnCount); // Compute the information to fill this vector. index_type lastIndex = 0; index_type currentNumberOfSetBits = 0; - - // If we are requested to add missing diagonal entries, we need to make sure the corresponding rows are also - // taken. -// storm::storage::BitVector columnBitCountConstraint = columnConstraint; -// if (insertDiagonalEntries) { -// columnBitCountConstraint |= rowGroupConstraint; -// } for (auto index : columnConstraint) { while (lastIndex <= index) { - bitsSetBeforeIndex.push_back(currentNumberOfSetBits); + columnBitsSetBeforeIndex.push_back(currentNumberOfSetBits); ++lastIndex; } ++currentNumberOfSetBits; } + std::vector* rowBitsSetBeforeIndex; + if (&rowGroupConstraint == &columnConstraint) { + rowBitsSetBeforeIndex = &columnBitsSetBeforeIndex; + } else { + rowBitsSetBeforeIndex = new std::vector(rowCount); + + lastIndex = 0; + currentNumberOfSetBits = 0; + for (auto index : rowGroupConstraint) { + while (lastIndex <= index) { + rowBitsSetBeforeIndex->push_back(currentNumberOfSetBits); + ++lastIndex; + } + ++currentNumberOfSetBits; + } + } // Then, we need to determine the number of entries and the number of rows of the submatrix. index_type subEntries = 0; @@ -551,7 +560,7 @@ namespace storm { if (columnConstraint.get(it->getColumn())) { ++subEntries; - if (index == it->getColumn()) { + if (columnBitsSetBeforeIndex[it->getColumn()] == (*rowBitsSetBeforeIndex)[index]) { foundDiagonalElement = true; } } @@ -580,13 +589,13 @@ namespace storm { for (const_iterator it = this->begin(i), ite = this->end(i); it != ite; ++it) { if (columnConstraint.get(it->getColumn())) { - if (index == it->getColumn()) { + if (columnBitsSetBeforeIndex[it->getColumn()] == (*rowBitsSetBeforeIndex)[index]) { insertedDiagonalElement = true; - } else if (insertDiagonalEntries && !insertedDiagonalElement && it->getColumn() > index) { - matrixBuilder.addNextValue(rowCount, rowCount, storm::utility::zero()); + } else if (insertDiagonalEntries && !insertedDiagonalElement && columnBitsSetBeforeIndex[it->getColumn()] > (*rowBitsSetBeforeIndex)[index]) { + matrixBuilder.addNextValue(rowCount, rowGroupCount, storm::utility::zero()); insertedDiagonalElement = true; } - matrixBuilder.addNextValue(rowCount, bitsSetBeforeIndex[it->getColumn()], it->getValue()); + matrixBuilder.addNextValue(rowCount, columnBitsSetBeforeIndex[it->getColumn()], it->getValue()); } } if (insertDiagonalEntries && !insertedDiagonalElement && rowGroupCount < submatrixColumnCount) { @@ -597,6 +606,11 @@ namespace storm { ++rowGroupCount; } + // If the constraints were not the same, we have allocated some additional memory we need to free now. + if (&rowGroupConstraint != &columnConstraint) { + delete rowBitsSetBeforeIndex; + } + return matrixBuilder.build(); } @@ -655,12 +669,12 @@ namespace storm { SparseMatrix SparseMatrix::transpose(bool joinGroups, bool keepZeros) const { index_type rowCount = this->getColumnCount(); index_type columnCount = joinGroups ? this->getRowGroupCount() : this->getRowCount(); - if (keepZeros) { - index_type entryCount = this->getEntryCount(); - } else { - this->updateNonzeroEntryCount(); - index_type entryCount = this->getNonzeroEntryCount(); - } + if (keepZeros) { + index_type entryCount = this->getEntryCount(); + } else { + this->updateNonzeroEntryCount(); + index_type entryCount = this->getNonzeroEntryCount(); + } std::vector rowIndications(rowCount + 1); std::vector> columnsAndValues(entryCount); @@ -668,7 +682,7 @@ namespace storm { // First, we need to count how many entries each column has. for (index_type group = 0; group < columnCount; ++group) { for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) { - if (transition.getValue() != storm::utility::zero() || keepZeros) { + if (transition.getValue() != storm::utility::zero() || keepZeros) { ++rowIndications[transition.getColumn() + 1]; } } @@ -687,7 +701,7 @@ namespace storm { // Now we are ready to actually fill in the values of the transposed matrix. for (index_type group = 0; group < columnCount; ++group) { for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) { - if (transition.getValue() != storm::utility::zero() || keepZeros) { + if (transition.getValue() != storm::utility::zero() || keepZeros) { columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue()); nextIndices[transition.getColumn()]++; } @@ -700,10 +714,10 @@ namespace storm { } storm::storage::SparseMatrix transposedMatrix(columnCount, std::move(rowIndications), std::move(columnsAndValues), std::move(rowGroupIndices), false); - + return transposedMatrix; } - + template void SparseMatrix::convertToEquationSystem() { invertDiagonal(); @@ -713,22 +727,22 @@ namespace storm { template void SparseMatrix::invertDiagonal() { // Now iterate over all row groups and set the diagonal elements to the inverted value. - // If there is a row without the diagonal element, an exception is thrown. - ValueType one = storm::utility::one(); - ValueType zero = storm::utility::zero(); + // If there is a row without the diagonal element, an exception is thrown. + ValueType one = storm::utility::one(); + ValueType zero = storm::utility::zero(); bool foundDiagonalElement = false; for (index_type group = 0; group < this->getRowGroupCount(); ++group) { for (auto& entry : this->getRowGroup(group)) { if (entry.getColumn() == group) { - if (entry.getValue() == one) { - --this->nonzeroEntryCount; - entry.setValue(zero); - } else if (entry.getValue() == zero) { - ++this->nonzeroEntryCount; - entry.setValue(one); - } else { - entry.setValue(one - entry.getValue()); - } + if (entry.getValue() == one) { + --this->nonzeroEntryCount; + entry.setValue(zero); + } else if (entry.getValue() == zero) { + ++this->nonzeroEntryCount; + entry.setValue(one); + } else { + entry.setValue(one - entry.getValue()); + } foundDiagonalElement = true; } } @@ -758,7 +772,7 @@ namespace storm { for (index_type group = 0; group < this->getRowGroupCount(); ++group) { for (auto& entry : this->getRowGroup(group)) { if (entry.getColumn() == group) { - --this->nonzeroEntryCount; + --this->nonzeroEntryCount; entry.setValue(storm::utility::zero()); } } @@ -839,15 +853,28 @@ namespace storm { typename std::vector::iterator resultIterator = result.begin(); typename std::vector::iterator resultIteratorEnd = result.end(); - for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { - *resultIterator = storm::utility::zero(); - - for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { - *resultIterator += it->getValue() * vector[it->getColumn()]; + // If the vector to multiply with and the target vector are actually the same, we need an auxiliary variable + // to store the intermediate result. + if (&vector == &result) { + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { + ValueType tmpValue = storm::utility::zero(); + + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { + tmpValue += it->getValue() * vector[it->getColumn()]; + } + *resultIterator = tmpValue; + } + } else { + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator) { + *resultIterator = storm::utility::zero(); + + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { + *resultIterator += it->getValue() * vector[it->getColumn()]; + } } } } - + #ifdef STORM_HAVE_INTELTBB template void SparseMatrix::multiplyWithVectorParallel(std::vector const& vector, std::vector& result) const { @@ -872,7 +899,37 @@ namespace storm { }); } #endif - + + template + void SparseMatrix::performSuccessiveOverRelaxationStep(ValueType omega, std::vector& x, std::vector const& b) const { + const_iterator it = this->begin(); + const_iterator ite; + std::vector::const_iterator rowIterator = rowIndications.begin(); + typename std::vector::const_iterator bIt = b.begin(); + typename std::vector::const_iterator bIte = b.end(); + typename std::vector::iterator resultIterator = x.begin(); + typename std::vector::iterator resultIteratorEnd = x.end(); + + // If the vector to multiply with and the target vector are actually the same, we need an auxiliary variable + // to store the intermediate result. + index_type currentRow = 0; + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++bIt) { + ValueType tmpValue = storm::utility::zero(); + ValueType diagonalElement = storm::utility::zero(); + + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { + if (it->getColumn() != currentRow) { + tmpValue += it->getValue() * x[it->getColumn()]; + } else { + diagonalElement += it->getValue(); + } + } + + *resultIterator = ((storm::utility::one() - omega) * *resultIterator) + (omega / diagonalElement) * (*bIt - tmpValue); + ++currentRow; + } + } + template void SparseMatrix::multiplyVectorWithMatrix(std::vector const& vector, std::vector& result) const { const_iterator it = this->begin(); @@ -982,7 +1039,7 @@ namespace storm { if (this->getRowCount() != matrix.getRowCount()) return false; if (this->getColumnCount() != matrix.getColumnCount()) return false; if (this->getRowGroupIndices() != matrix.getRowGroupIndices()) return false; - + // Check the subset property for all rows individually. for (index_type row = 0; row < this->getRowCount(); ++row) { for (const_iterator it1 = this->begin(row), ite1 = this->end(row), it2 = matrix.begin(row), ite2 = matrix.end(row); it1 != ite1; ++it1) { @@ -1054,21 +1111,21 @@ namespace storm { } // Explicitly instantiate the entry, builder and the matrix. - //double - template class MatrixEntry::index_type, double>; - template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); - template class SparseMatrixBuilder; - template class SparseMatrix; - template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); - //float - template class MatrixEntry::index_type, float>; - template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); - template class SparseMatrixBuilder; - template class SparseMatrix; - template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); - //int - template class MatrixEntry::index_type, int>; - template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); + //double + template class MatrixEntry::index_type, double>; + template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); + template class SparseMatrixBuilder; + template class SparseMatrix; + template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); + //float + template class MatrixEntry::index_type, float>; + template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); + template class SparseMatrixBuilder; + template class SparseMatrix; + template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); + //int + template class MatrixEntry::index_type, int>; + template std::ostream& operator<<(std::ostream& out, MatrixEntry const& entry); template class SparseMatrixBuilder; template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 33863da48..df682d756 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -692,6 +692,15 @@ namespace storm { */ void multiplyVectorWithMatrix(std::vector const& vector, std::vector& result) const; + /*! + * Performs one step of the successive over-relaxation technique. + * + * @param omega The Omega parameter for SOR. + * @param x The current solution vector. The result will be written to the very same vector. + * @param b The 'right-hand side' of the problem. + */ + void performSuccessiveOverRelaxationStep(ValueType omega, std::vector& x, std::vector const& b) const; + /*! * Multiplies the matrix with the given vector in a sequential way and writes the result to the given result * vector. diff --git a/src/utility/solver.cpp b/src/utility/solver.cpp index 2be080beb..785f1d00c 100644 --- a/src/utility/solver.cpp +++ b/src/utility/solver.cpp @@ -41,9 +41,28 @@ namespace storm { return std::unique_ptr>(new storm::solver::GmmxxLinearEquationSolver(matrix)); } + template + NativeLinearEquationSolverFactory::NativeLinearEquationSolverFactory() { + switch (storm::settings::nativeEquationSolverSettings().getLinearEquationSystemMethod()) { + case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Jacobi: + this->method = storm::solver::NativeLinearEquationSolver::SolutionMethod::Jacobi; + break; + case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel: + this->method = storm::solver::NativeLinearEquationSolver::SolutionMethod::GaussSeidel; + case settings::modules::NativeEquationSolverSettings::LinearEquationMethod::SOR: + this->method = storm::solver::NativeLinearEquationSolver::SolutionMethod::SOR; + } + omega = storm::settings::nativeEquationSolverSettings().getOmega(); + } + + template + NativeLinearEquationSolverFactory::NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolver::SolutionMethod method, ValueType omega) : method(method), omega(omega) { + // Intentionally left empty. + } + template std::unique_ptr> NativeLinearEquationSolverFactory::create(storm::storage::SparseMatrix const& matrix) const { - return std::unique_ptr>(new storm::solver::NativeLinearEquationSolver(matrix)); + return std::unique_ptr>(new storm::solver::NativeLinearEquationSolver(matrix, method)); } template diff --git a/src/utility/solver.h b/src/utility/solver.h index 38cc7fbc6..cd56695f0 100644 --- a/src/utility/solver.h +++ b/src/utility/solver.h @@ -5,10 +5,12 @@ #include "src/solver/SymbolicLinearEquationSolver.h" #include "src/solver/LinearEquationSolver.h" +#include "src/solver/NativeLinearEquationSolver.h" #include "src/solver/MinMaxLinearEquationSolver.h" #include "src/solver/LpSolver.h" #include "src/storage/dd/DdType.h" +#include "src/settings/modules/NativeEquationSolverSettings.h" #include "src/exceptions/InvalidSettingsException.h" @@ -42,7 +44,14 @@ namespace storm { template class NativeLinearEquationSolverFactory : public LinearEquationSolverFactory { public: + NativeLinearEquationSolverFactory(); + NativeLinearEquationSolverFactory(typename storm::solver::NativeLinearEquationSolver::SolutionMethod method, ValueType omega); + virtual std::unique_ptr> create(storm::storage::SparseMatrix const& matrix) const override; + + private: + typename storm::solver::NativeLinearEquationSolver::SolutionMethod method; + ValueType omega; }; template diff --git a/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp index 387f54ad4..f58b4259e 100644 --- a/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp +++ b/test/functional/modelchecker/NativeDtmcPrctlModelCheckerTest.cpp @@ -1,6 +1,7 @@ #include "gtest/gtest.h" #include "storm-config.h" +#include "src/settings/SettingMemento.h" #include "src/logic/Formulas.h" #include "src/utility/solver.h" #include "src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h" @@ -144,7 +145,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) { dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory(storm::solver::NativeLinearEquationSolver::SolutionMethod::SOR, 0.9))); auto labelFormula = std::make_shared("a"); auto lraFormula = std::make_shared(labelFormula); @@ -169,7 +170,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) { dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory(storm::solver::NativeLinearEquationSolver::SolutionMethod::SOR, 0.9))); auto labelFormula = std::make_shared("a"); auto lraFormula = std::make_shared(labelFormula); @@ -194,7 +195,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRASingleBscc) { dtmc.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + storm::modelchecker::SparseDtmcPrctlModelChecker checker(*dtmc, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory(storm::solver::NativeLinearEquationSolver::SolutionMethod::SOR, 0.9))); auto labelFormula = std::make_shared("a"); auto lraFormula = std::make_shared(labelFormula); @@ -255,7 +256,7 @@ TEST(SparseDtmcPrctlModelCheckerTest, LRA) { mdp.reset(new storm::models::sparse::Dtmc(transitionMatrix, ap, boost::none, boost::none, boost::none)); - storm::modelchecker::SparseDtmcPrctlModelChecker checker(*mdp, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory())); + storm::modelchecker::SparseDtmcPrctlModelChecker checker(*mdp, std::unique_ptr>(new storm::utility::solver::NativeLinearEquationSolverFactory(storm::solver::NativeLinearEquationSolver::SolutionMethod::SOR, 0.9))); auto labelFormula = std::make_shared("a"); auto lraFormula = std::make_shared(labelFormula);