Browse Source

Merge branch 'future' into multi-objective

Former-commit-id: e59a5483cc
tempestpy_adaptions
TimQu 9 years ago
parent
commit
98898dde84
  1. 70
      src/builder/ExplicitModelBuilder.cpp
  2. 19
      src/builder/ExplicitModelBuilder.h
  3. 4
      src/generator/PrismNextStateGenerator.cpp
  4. 9
      src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
  5. 7
      src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
  6. 6
      src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp
  7. 2
      src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp
  8. 12
      src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp
  9. 51
      src/models/sparse/MarkovAutomaton.cpp
  10. 41
      src/models/sparse/MarkovAutomaton.h
  11. 31
      src/solver/EigenLinearEquationSolver.cpp
  12. 12
      src/solver/EigenLinearEquationSolver.h
  13. 42
      src/solver/EliminationLinearEquationSolver.cpp
  14. 16
      src/solver/EliminationLinearEquationSolver.h
  15. 121
      src/solver/GmmxxLinearEquationSolver.cpp
  16. 26
      src/solver/GmmxxLinearEquationSolver.h
  17. 72
      src/solver/LinearEquationSolver.cpp
  18. 72
      src/solver/LinearEquationSolver.h
  19. 23
      src/solver/MinMaxLinearEquationSolver.cpp
  20. 45
      src/solver/MinMaxLinearEquationSolver.h
  21. 153
      src/solver/NativeLinearEquationSolver.cpp
  22. 24
      src/solver/NativeLinearEquationSolver.h
  23. 172
      src/solver/StandardMinMaxLinearEquationSolver.cpp
  24. 19
      src/solver/StandardMinMaxLinearEquationSolver.h
  25. 18
      src/solver/TopologicalMinMaxLinearEquationSolver.cpp
  26. 4
      src/solver/TopologicalMinMaxLinearEquationSolver.h
  27. 10
      src/storage/BitVector.cpp
  28. 11
      src/storage/BitVector.h
  29. 64
      src/storage/SparseMatrix.cpp
  30. 8
      src/storage/SparseMatrix.h
  31. 2
      src/utility/ModelInstantiator.h
  32. 26
      test/functional/builder/ExplicitPrismModelBuilderTest.cpp
  33. 19
      test/functional/builder/hybrid_states.ma
  34. 15
      test/functional/builder/simple.ma
  35. 50
      test/functional/builder/stream2.ma
  36. 10
      test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp
  37. 10
      test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp

70
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<storm::models::sparse::Model<ValueType, RewardModelType>>(new storm::models::sparse::Mdp<ValueType, RewardModelType>(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<storm::models::sparse::Model<ValueType, RewardModelType>>(new storm::models::sparse::MarkovAutomaton<ValueType, RewardModelType>(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 <typename ValueType, typename RewardModelType, typename StateType>
boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> ExplicitModelBuilder<ValueType, RewardModelType, StateType>::buildMatrices(storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<RewardModelBuilder<typename RewardModelType::ValueType>>& rewardModelBuilders) {
void ExplicitModelBuilder<ValueType, RewardModelType, StateType>::buildMatrices(storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<RewardModelBuilder<typename RewardModelType::ValueType>>& rewardModelBuilders, boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>& choiceLabels , boost::optional<storm::storage::BitVector>& markovianChoices) {
// Create choice labels, if requested,
boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> choiceLabels;
if (generator->getOptions().isBuildChoiceLabelsSet()) {
choiceLabels = std::vector<boost::container::flat_set<uint_fast64_t>>();
}
// 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<StateType (CompressedState const&)> stateToIdCallback = std::bind(&ExplicitModelBuilder<ValueType, RewardModelType, StateType>::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<uint_fast64_t>());
}
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 <typename ValueType, typename RewardModelType, typename StateType>
@ -323,7 +347,8 @@ namespace storm {
rewardModelBuilders.emplace_back(generator->getRewardModelInformation(i));
}
modelComponents.choiceLabeling = buildMatrices(transitionMatrixBuilder, rewardModelBuilders);
boost::optional<storm::storage::BitVector> 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 <typename ValueType, typename RewardModelType, typename StateType>
void ExplicitModelBuilder<ValueType, RewardModelType, StateType>::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 <typename ValueType, typename RewardModelType, typename StateType>
storm::models::sparse::StateLabeling ExplicitModelBuilder<ValueType, RewardModelType, StateType>::buildStateLabeling() {
return generator->label(stateStorage.stateToId, stateStorage.initialStateIndices, stateStorage.deadlockStateIndices);

19
src/builder/ExplicitModelBuilder.h

@ -63,6 +63,9 @@ namespace storm {
// A vector that stores a labeling for each choice.
boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> choiceLabeling;
// A vector that stores which states are markovian.
boost::optional<storm::storage::BitVector> 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<std::vector<boost::container::flat_set<uint_fast64_t>>> buildMatrices(storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<RewardModelBuilder<typename RewardModelType::ValueType>>& rewardModelBuilders);
void buildMatrices(storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<RewardModelBuilder<typename RewardModelType::ValueType>>& rewardModelBuilders, boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>>& choiceLabels , boost::optional<storm::storage::BitVector>& 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.
*

4
src/generator/PrismNextStateGenerator.cpp

@ -366,7 +366,7 @@ namespace storm {
continue;
}
result.push_back(Choice<ValueType>(command.getActionIndex()));
result.push_back(Choice<ValueType>(command.getActionIndex(), command.isMarkovian()));
Choice<ValueType>& choice = result.back();
// Remember the command labels only if we were asked to.
@ -564,4 +564,4 @@ namespace storm {
template class PrismNextStateGenerator<storm::RationalNumber>;
template class PrismNextStateGenerator<storm::RationalFunction>;
}
}
}

9
src/modelchecker/csl/helper/SparseCtmcCslHelper.cpp

@ -627,19 +627,19 @@ namespace storm {
result = std::vector<ValueType>(values.size());
}
}
std::vector<ValueType> multiplicationResult(result.size());
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(uniformizedMatrix));
solver->allocateAuxMemory(storm::solver::LinearEquationSolverOperation::MultiplyRepeatedly);
if (!useMixedPoissonProbabilities && std::get<0>(foxGlynnResult) > 1) {
// Perform the matrix-vector multiplications (without adding).
solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1, &multiplicationResult);
solver->repeatedMultiply(values, addVector, std::get<0>(foxGlynnResult) - 1);
} else if (useMixedPoissonProbabilities) {
std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&uniformizationRate] (ValueType const& a, ValueType const& b) { return a + b / uniformizationRate; };
// For the iterations below the left truncation point, we need to add and scale the result with the uniformization rate.
for (uint_fast64_t index = 1; index < startingIteration; ++index) {
solver->repeatedMultiply(values, nullptr, 1, &multiplicationResult);
solver->repeatedMultiply(values, nullptr, 1);
storm::utility::vector::applyPointwise(result, values, result, addAndScale);
}
}
@ -649,7 +649,7 @@ namespace storm {
ValueType weight = 0;
std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; };
for (uint_fast64_t index = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) {
solver->repeatedMultiply(values, addVector, 1, &multiplicationResult);
solver->repeatedMultiply(values, addVector, 1);
weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
storm::utility::vector::applyPointwise(result, values, result, addAndScale);
@ -658,7 +658,6 @@ namespace storm {
return result;
}
template <typename ValueType>
storm::storage::SparseMatrix<ValueType> SparseCtmcCslHelper::computeProbabilityMatrix(storm::storage::SparseMatrix<ValueType> const& rateMatrix, std::vector<ValueType> const& exitRates) {
// Turn the rates into probabilities by scaling each row with the exit rate of the state.

7
src/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp

@ -86,6 +86,7 @@ namespace storm {
}
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(aProbabilistic);
solver->allocateAuxMemory(storm::solver::MinMaxLinearEquationSolverOperation::SolveEquations);
// Perform the actual value iteration
// * loop until the step bound has been reached
@ -94,15 +95,13 @@ namespace storm {
// and 1_G being the characteristic vector for all goal states.
// * perform one timed-step using v_MS := A_MSwG * v_MS + A_MStoPS * v_PS + (A * 1_G)|MS
std::vector<ValueType> markovianNonGoalValuesSwap(markovianNonGoalValues);
std::vector<ValueType> multiplicationResultScratchMemory(aProbabilistic.getRowCount());
std::vector<ValueType> aProbabilisticScratchMemory(probabilisticNonGoalValues.size());
for (uint_fast64_t currentStep = 0; currentStep < numberOfSteps; ++currentStep) {
// Start by (re-)computing bProbabilistic = bProbabilisticFixed + aProbabilisticToMarkovian * vMarkovian.
aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
// Now perform the inner value iteration for probabilistic states.
solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic);
// (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic.
aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian);
@ -115,7 +114,7 @@ namespace storm {
// After the loop, perform one more step of the value iteration for PS states.
aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
solver->solveEquations(dir, probabilisticNonGoalValues, bProbabilistic);
}
template<typename ValueType>

6
src/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp

@ -141,7 +141,7 @@ namespace storm {
std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd);
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first));
solver->multiply(dir, x, &explicitRepresentation.second, stepBound);
solver->repeatedMultiply(dir, x, &explicitRepresentation.second, stepBound);
// Return a hybrid check result that stores the numerical values explicitly.
return std::unique_ptr<CheckResult>(new storm::modelchecker::HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getReachableStates() && !maybeStates, psiStates.template toAdd<ValueType>(), maybeStates, odd, x));
@ -166,7 +166,7 @@ namespace storm {
// Perform the matrix-vector multiplication.
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitMatrix));
solver->multiply(dir, x, nullptr, stepBound);
solver->repeatedMultiply(dir, x, nullptr, stepBound);
// Return a hybrid check result that stores the numerical values explicitly.
return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero<ValueType>(), model.getReachableStates(), odd, x));
@ -195,7 +195,7 @@ namespace storm {
// Perform the matrix-vector multiplication.
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first));
solver->multiply(dir, x, &explicitRepresentation.second, stepBound);
solver->repeatedMultiply(dir, x, &explicitRepresentation.second, stepBound);
// Return a hybrid check result that stores the numerical values explicitly.
return std::unique_ptr<CheckResult>(new HybridQuantitativeCheckResult<DdType>(model.getReachableStates(), model.getManager().getBddZero(), model.getManager().template getAddZero<ValueType>(), model.getReachableStates(), odd, x));

2
src/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp

@ -121,7 +121,7 @@ namespace storm {
// Perform one single matrix-vector multiplication.
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(transitionMatrix);
solver->repeatedMultiply(result);
solver->repeatedMultiply(result, nullptr, 1);
return result;
}

12
src/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp

@ -49,7 +49,7 @@ namespace storm {
std::vector<ValueType> subresult(maybeStates.getNumberOfSetBits());
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(std::move(submatrix));
solver->multiply(dir, subresult, &b, stepBound);
solver->repeatedMultiply(dir, subresult, &b, stepBound);
// Set the values of the resulting vector accordingly.
storm::utility::vector::setVectorValues(result, maybeStates, subresult);
@ -67,7 +67,7 @@ namespace storm {
storm::utility::vector::setVectorValues(result, nextStates, storm::utility::one<ValueType>());
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
solver->multiply(dir, result);
solver->repeatedMultiply(dir, result, nullptr, 1);
return result;
}
@ -211,7 +211,7 @@ namespace storm {
std::vector<ValueType> result(rewardModel.getStateRewardVector());
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
solver->multiply(dir, result, nullptr, stepCount);
solver->repeatedMultiply(dir, result, nullptr, stepCount);
return result;
}
@ -235,7 +235,7 @@ namespace storm {
}
std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver = minMaxLinearEquationSolverFactory.create(transitionMatrix);
solver->multiply(dir, result, &totalRewardVector, stepBound);
solver->repeatedMultiply(dir, result, &totalRewardVector, stepBound);
return result;
}
@ -601,13 +601,13 @@ namespace storm {
if (fixedTargetStates.get(state)) {
builder.addNextValue(currentRow, newGoalState, conditionProbabilities[state]);
if (!storm::utility::isZero(conditionProbabilities[state])) {
builder.addNextValue(currentRow, newFailState, 1 - conditionProbabilities[state]);
builder.addNextValue(currentRow, newFailState, storm::utility::one<ValueType>() - conditionProbabilities[state]);
}
++currentRow;
} else if (conditionStates.get(state)) {
builder.addNextValue(currentRow, newGoalState, targetProbabilities[state]);
if (!storm::utility::isZero(targetProbabilities[state])) {
builder.addNextValue(currentRow, newStopState, 1 - targetProbabilities[state]);
builder.addNextValue(currentRow, newStopState, storm::utility::one<ValueType>() - targetProbabilities[state]);
}
++currentRow;
} else {

51
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<ValueType, RewardModelType>(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 <typename ValueType, typename RewardModelType>
MarkovAutomaton<ValueType, RewardModelType>::MarkovAutomaton(storm::storage::SparseMatrix<ValueType> const& rateMatrix,
storm::models::sparse::StateLabeling const& stateLabeling,
storm::storage::BitVector const& markovianStates,
std::unordered_map<std::string, RewardModelType> const& rewardModels,
boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling)
: NondeterministicModel<ValueType, RewardModelType>(storm::models::ModelType::MarkovAutomaton, rateMatrix, stateLabeling, rewardModels, optionalChoiceLabeling), markovianStates(markovianStates) {
turnRatesToProbabilities();
this->closed = checkIsClosed();
}
template <typename ValueType, typename RewardModelType>
MarkovAutomaton<ValueType, RewardModelType>::MarkovAutomaton(storm::storage::SparseMatrix<ValueType>&& rateMatrix,
storm::models::sparse::StateLabeling&& stateLabeling,
storm::storage::BitVector&& markovianStates,
std::unordered_map<std::string, RewardModelType>&& rewardModels,
boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling)
: NondeterministicModel<ValueType, RewardModelType>(storm::models::ModelType::MarkovAutomaton, rateMatrix, stateLabeling, rewardModels, optionalChoiceLabeling), markovianStates(markovianStates) {
turnRatesToProbabilities();
this->closed = checkIsClosed();
}
template <typename ValueType, typename RewardModelType>
MarkovAutomaton<ValueType, RewardModelType>::MarkovAutomaton(storm::storage::SparseMatrix<ValueType>&& transitionMatrix,
@ -230,9 +256,18 @@ namespace storm {
template <typename ValueType, typename RewardModelType>
void MarkovAutomaton<ValueType, RewardModelType>::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.");
}
}
}

41
src/models/sparse/MarkovAutomaton.h

@ -31,7 +31,43 @@ namespace storm {
std::vector<ValueType> const& exitRates,
std::unordered_map<std::string, RewardModelType> const& rewardModels = std::unordered_map<std::string, RewardModelType>(),
boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
/*!
* 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<ValueType> const& rateMatrix,
storm::models::sparse::StateLabeling const& stateLabeling,
storm::storage::BitVector const& markovianStates,
std::unordered_map<std::string, RewardModelType> const& rewardModels = std::unordered_map<std::string, RewardModelType>(),
boost::optional<std::vector<LabelSet>> const& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
/*!
* 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<ValueType>&& rateMatrix,
storm::models::sparse::StateLabeling&& stateLabeling,
storm::storage::BitVector&& markovianStates,
std::unordered_map<std::string, RewardModelType>&& rewardModels = std::unordered_map<std::string, RewardModelType>(),
boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
/*!
* Constructs a model by moving the given data.
*
@ -48,7 +84,7 @@ namespace storm {
std::vector<ValueType> const& exitRates,
std::unordered_map<std::string, RewardModelType>&& rewardModels = std::unordered_map<std::string, RewardModelType>(),
boost::optional<std::vector<LabelSet>>&& optionalChoiceLabeling = boost::optional<std::vector<LabelSet>>());
/*!
* 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();

31
src/solver/EigenLinearEquationSolver.cpp

@ -109,12 +109,23 @@ namespace storm {
template<typename ValueType>
EigenLinearEquationSolver<ValueType>::EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EigenLinearEquationSolverSettings<ValueType> const& settings) : settings(settings) {
this->setMatrix(std::move(A));
}
template<typename ValueType>
void EigenLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(A);
}
template<typename ValueType>
void EigenLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
// Take ownership of the matrix so it is destroyed after we have translated it to Eigen's format.
storm::storage::SparseMatrix<ValueType> localA(std::move(A));
eigenA = storm::adapters::EigenAdapter::toEigenSparseMatrix<ValueType>(localA);
this->setMatrix(localA);
}
template<typename ValueType>
void EigenLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
void EigenLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// Map the input vectors to Eigen's format.
auto eigenX = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b.data(), b.size());
@ -197,7 +208,7 @@ namespace storm {
}
template<typename ValueType>
void EigenLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
void EigenLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
// Typedef the map-type so we don't have to spell it out.
typedef decltype(Eigen::Matrix<ValueType, Eigen::Dynamic, 1>::Map(b->data(), b->size())) MapType;
@ -234,10 +245,20 @@ namespace storm {
return settings;
}
template<typename ValueType>
uint64_t EigenLinearEquationSolver<ValueType>::getMatrixRowCount() const {
return eigenA->rows();
}
template<typename ValueType>
uint64_t EigenLinearEquationSolver<ValueType>::getMatrixColumnCount() const {
return eigenA->cols();
}
// Specialization form storm::RationalNumber
template<>
void EigenLinearEquationSolver<storm::RationalNumber>::solveEquations(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b, std::vector<storm::RationalNumber>* multiplyResult) const {
void EigenLinearEquationSolver<storm::RationalNumber>::solveEquations(std::vector<storm::RationalNumber>& x, std::vector<storm::RationalNumber> const& b) const {
// Map the input vectors to Eigen's format.
auto eigenX = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<storm::RationalNumber, Eigen::Dynamic, 1>::Map(b.data(), b.size());
@ -250,7 +271,7 @@ namespace storm {
// Specialization form storm::RationalFunction
template<>
void EigenLinearEquationSolver<storm::RationalFunction>::solveEquations(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b, std::vector<storm::RationalFunction>* multiplyResult) const {
void EigenLinearEquationSolver<storm::RationalFunction>::solveEquations(std::vector<storm::RationalFunction>& x, std::vector<storm::RationalFunction> const& b) const {
// Map the input vectors to Eigen's format.
auto eigenX = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(x.data(), x.size());
auto eigenB = Eigen::Matrix<storm::RationalFunction, Eigen::Dynamic, 1>::Map(b.data(), b.size());

12
src/solver/EigenLinearEquationSolver.h

@ -61,13 +61,19 @@ namespace storm {
EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EigenLinearEquationSolverSettings<ValueType> const& settings = EigenLinearEquationSolverSettings<ValueType>());
EigenLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EigenLinearEquationSolverSettings<ValueType> const& settings = EigenLinearEquationSolverSettings<ValueType>());
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
EigenLinearEquationSolverSettings<ValueType>& getSettings();
EigenLinearEquationSolverSettings<ValueType> const& getSettings() const;
private:
virtual uint64_t getMatrixRowCount() const override;
virtual uint64_t getMatrixColumnCount() const override;
// The (eigen) matrix associated with this equation solver.
std::unique_ptr<Eigen::SparseMatrix<ValueType>> eigenA;

42
src/solver/EliminationLinearEquationSolver.cpp

@ -34,19 +34,29 @@ namespace storm {
}
template<typename ValueType>
EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), settings(settings) {
// Intentionally left empty.
EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings) {
this->setMatrix(A);
}
template<typename ValueType>
EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), settings(settings) {
// Intentionally left empty.
EliminationLinearEquationSolver<ValueType>::EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings) {
this->setMatrix(std::move(A));
}
template<typename ValueType>
void EliminationLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
STORM_LOG_WARN_COND(multiplyResult == nullptr, "Providing scratch memory will not improve the performance of this solver.");
void EliminationLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
this->A = &A;
localA.reset();
}
template<typename ValueType>
void EliminationLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
localA = std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A));
this->A = localA.get();
}
template<typename ValueType>
void EliminationLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// FIXME: This solver will not work for all input systems. More concretely, the current implementation will
// not work for systems that have a 0 on the diagonal. This is not a restriction of this technique in general
// but arbitrary matrices require pivoting, which is not currently implemented.
@ -59,7 +69,7 @@ namespace storm {
if (localA) {
localA->convertToEquationSystem();
} else {
locallyConvertedMatrix = A;
locallyConvertedMatrix = *A;
locallyConvertedMatrix.convertToEquationSystem();
}
@ -101,16 +111,16 @@ namespace storm {
}
template<typename ValueType>
void EliminationLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
void EliminationLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
if (&x != &result) {
A.multiplyWithVector(x, result);
A->multiplyWithVector(x, result);
if (b != nullptr) {
storm::utility::vector::addVectors(result, *b, result);
}
} else {
// If the two vectors are aliases, we need to create a temporary.
std::vector<ValueType> tmp(result.size());
A.multiplyWithVector(x, tmp);
A->multiplyWithVector(x, tmp);
if (b != nullptr) {
storm::utility::vector::addVectors(tmp, *b, result);
}
@ -127,6 +137,16 @@ namespace storm {
return settings;
}
template<typename ValueType>
uint64_t EliminationLinearEquationSolver<ValueType>::getMatrixRowCount() const {
return this->A->getRowCount();
}
template<typename ValueType>
uint64_t EliminationLinearEquationSolver<ValueType>::getMatrixColumnCount() const {
return this->A->getColumnCount();
}
template<typename ValueType>
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> EliminationLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
return std::make_unique<storm::solver::EliminationLinearEquationSolver<ValueType>>(matrix, settings);

16
src/solver/EliminationLinearEquationSolver.h

@ -29,8 +29,11 @@ namespace storm {
EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, EliminationLinearEquationSolverSettings<ValueType> const& settings = EliminationLinearEquationSolverSettings<ValueType>());
EliminationLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, EliminationLinearEquationSolverSettings<ValueType> const& settings = EliminationLinearEquationSolverSettings<ValueType>());
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
EliminationLinearEquationSolverSettings<ValueType>& getSettings();
EliminationLinearEquationSolverSettings<ValueType> const& getSettings() const;
@ -38,13 +41,16 @@ namespace storm {
private:
void initializeSettings();
virtual uint64_t getMatrixRowCount() const override;
virtual uint64_t getMatrixColumnCount() const override;
// If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
// when the solver is destructed.
std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
// A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
// the reference refers to localA.
storm::storage::SparseMatrix<ValueType> const& A;
// A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix
// the pointer refers to localA.
storm::storage::SparseMatrix<ValueType> const* A;
// The settings used by the solver.
EliminationLinearEquationSolverSettings<ValueType> settings;

121
src/solver/GmmxxLinearEquationSolver.cpp

@ -111,17 +111,31 @@ namespace storm {
}
template<typename ValueType>
GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A)), settings(settings) {
// Intentionally left empty.
GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiMemory(nullptr) {
this->setMatrix(A);
}
template<typename ValueType>
GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), gmmxxMatrix(storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(*localA)), settings(settings) {
// Intentionally left empty.
GmmxxLinearEquationSolver<ValueType>::GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), gmmxxMatrix(nullptr), settings(settings), auxiliaryJacobiMemory(nullptr) {
this->setMatrix(std::move(A));
}
template<typename ValueType>
void GmmxxLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
void GmmxxLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
localA.reset();
this->A = &A;
gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(A);
}
template<typename ValueType>
void GmmxxLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
localA = std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A));
this->A = localA.get();
gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(*localA);
}
template<typename ValueType>
void GmmxxLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
auto method = this->getSettings().getSolutionMethod();
auto preconditioner = this->getSettings().getPreconditioner();
STORM_LOG_INFO("Using method '" << method << "' with preconditioner '" << preconditioner << "' (max. " << this->getSettings().getMaximalNumberOfIterations() << " iterations).");
@ -166,7 +180,7 @@ namespace storm {
STORM_LOG_WARN("Iterative solver did not converge.");
}
} else if (method == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(A, x, b, multiplyResult);
uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(*A, x, b);
// Check if the solver converged and issue a warning otherwise.
if (iterations < this->getSettings().getMaximalNumberOfIterations()) {
@ -178,7 +192,7 @@ namespace storm {
}
template<typename ValueType>
void GmmxxLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
void GmmxxLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
if (b) {
gmm::mult_add(*gmmxxMatrix, x, *b, result);
} else {
@ -187,25 +201,20 @@ namespace storm {
}
template<typename ValueType>
uint_fast64_t GmmxxLinearEquationSolver<ValueType>::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
uint_fast64_t GmmxxLinearEquationSolver<ValueType>::solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
bool allocatedAuxMemory = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations);
if (allocatedAuxMemory) {
this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
}
// Get a Jacobi decomposition of the matrix A.
std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
// Convert the LU matrix to gmm++'s format.
std::unique_ptr<gmm::csr_matrix<ValueType>> gmmLU = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.first));
// To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
bool multiplyResultProvided = true;
std::vector<ValueType>* nextX = multiplyResult;
if (nextX == nullptr) {
nextX = new std::vector<ValueType>(x.size());
multiplyResultProvided = false;
}
std::vector<ValueType> const* copyX = nextX;
std::vector<ValueType>* currentX = &x;
// Target vector for precision calculation.
std::vector<ValueType> tmpX(x.size());
std::vector<ValueType>* nextX = auxiliaryJacobiMemory.get();
// Set up additional environment variables.
uint_fast64_t iterationCount = 0;
@ -213,9 +222,8 @@ namespace storm {
while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
// Compute D^-1 * (b - LU * x) and store result in nextX.
gmm::mult(*gmmLU, *currentX, tmpX);
gmm::add(b, gmm::scaled(tmpX, -storm::utility::one<ValueType>()), tmpX);
storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
gmm::mult_add(*gmmLU, gmm::scaled(*currentX, -storm::utility::one<ValueType>()), b, *nextX);
storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX);
// Now check if the process already converged within our precision.
converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion());
@ -229,13 +237,13 @@ namespace storm {
// If the last iteration did not write to the original x we have to swap the contents, because the
// output has to be written to the input parameter x.
if (currentX == copyX) {
if (currentX == auxiliaryJacobiMemory.get()) {
std::swap(x, *currentX);
}
// If the vector for the temporary multiplication result was not provided, we need to delete it.
if (!multiplyResultProvided) {
delete copyX;
// If we allocated auxiliary memory, we need to dispose of it now.
if (allocatedAuxMemory) {
this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
}
return iterationCount;
@ -251,6 +259,67 @@ namespace storm {
return settings;
}
template<typename ValueType>
bool GmmxxLinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::SolveEquations) {
if (this->getSettings().getSolutionMethod() == GmmxxLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
if (!auxiliaryJacobiMemory) {
auxiliaryJacobiMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount());
result = true;
}
}
}
result |= LinearEquationSolver<ValueType>::allocateAuxMemory(operation);
return result;
}
template<typename ValueType>
bool GmmxxLinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::SolveEquations) {
if (auxiliaryJacobiMemory) {
result = true;
auxiliaryJacobiMemory.reset();
}
}
result |= LinearEquationSolver<ValueType>::deallocateAuxMemory(operation);
return result;
}
template<typename ValueType>
bool GmmxxLinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::SolveEquations) {
if (auxiliaryJacobiMemory) {
result = auxiliaryJacobiMemory->size() != this->getMatrixColumnCount();
auxiliaryJacobiMemory->resize(this->getMatrixRowCount());
}
}
result |= LinearEquationSolver<ValueType>::reallocateAuxMemory(operation);
return result;
}
template<typename ValueType>
bool GmmxxLinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::SolveEquations) {
result |= static_cast<bool>(auxiliaryJacobiMemory);
}
result |= LinearEquationSolver<ValueType>::hasAuxMemory(operation);
return result;
}
template<typename ValueType>
uint64_t GmmxxLinearEquationSolver<ValueType>::getMatrixRowCount() const {
return this->A->getRowCount();
}
template<typename ValueType>
uint64_t GmmxxLinearEquationSolver<ValueType>::getMatrixColumnCount() const {
return this->A->getColumnCount();
}
template<typename ValueType>
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> GmmxxLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
return std::make_unique<storm::solver::GmmxxLinearEquationSolver<ValueType>>(matrix, settings);

26
src/solver/GmmxxLinearEquationSolver.h

@ -89,12 +89,20 @@ namespace storm {
GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings = GmmxxLinearEquationSolverSettings<ValueType>());
GmmxxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, GmmxxLinearEquationSolverSettings<ValueType> const& settings = GmmxxLinearEquationSolverSettings<ValueType>());
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
GmmxxLinearEquationSolverSettings<ValueType>& getSettings();
GmmxxLinearEquationSolverSettings<ValueType> const& getSettings() const;
virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const override;
virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const override;
virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const override;
virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const override;
private:
/*!
* Solves the linear equation system A*x = b given by the parameters using the Jacobi method.
@ -108,21 +116,27 @@ namespace storm {
* @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
* vector must be equal to the number of rows of A.
*/
uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const;
uint_fast64_t solveLinearEquationSystemWithJacobi(storm::storage::SparseMatrix<ValueType> const& A, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
virtual uint64_t getMatrixRowCount() const override;
virtual uint64_t getMatrixColumnCount() const override;
// If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
// when the solver is destructed.
std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
// A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
// the reference refers to localA.
storm::storage::SparseMatrix<ValueType> const& A;
// A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix
// the pointer refers to localA.
storm::storage::SparseMatrix<ValueType> const* A;
// The (gmm++) matrix associated with this equation solver.
std::unique_ptr<gmm::csr_matrix<ValueType>> gmmxxMatrix;
// The settings used by the solver.
GmmxxLinearEquationSolverSettings<ValueType> settings;
// Auxiliary storage for the Jacobi method.
mutable std::unique_ptr<std::vector<ValueType>> auxiliaryJacobiMemory;
};
template<typename ValueType>

72
src/solver/LinearEquationSolver.cpp

@ -14,36 +14,78 @@ namespace storm {
namespace solver {
template<typename ValueType>
void LinearEquationSolver<ValueType>::repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
LinearEquationSolver<ValueType>::LinearEquationSolver() : auxiliaryRepeatedMultiplyMemory(nullptr) {
// Intentionally left empty.
}
template<typename ValueType>
void LinearEquationSolver<ValueType>::repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n) const {
bool allocatedAuxMemory = !this->hasAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly);
if (allocatedAuxMemory) {
this->allocateAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly);
}
// Set up some temporary variables so that we can just swap pointers instead of copying the result after
// each iteration.
std::vector<ValueType>* currentX = &x;
bool multiplyResultProvided = true;
std::vector<ValueType>* nextX = multiplyResult;
if (nextX == nullptr) {
nextX = new std::vector<ValueType>(x.size());
multiplyResultProvided = false;
}
std::vector<ValueType> const* copyX = nextX;
std::vector<ValueType>* nextX = auxiliaryRepeatedMultiplyMemory.get();
// Now perform matrix-vector multiplication as long as we meet the bound.
for (uint_fast64_t i = 0; i < n; ++i) {
this->multiply(*currentX, *nextX, b);
this->multiply(*currentX, b, *nextX);
std::swap(nextX, currentX);
}
// If we performed an odd number of repetitions, we need to swap the contents of currentVector and x,
// because the output is supposed to be stored in the input vector x.
if (currentX == copyX) {
if (currentX == auxiliaryRepeatedMultiplyMemory.get()) {
std::swap(x, *currentX);
}
// If the vector for the temporary multiplication result was not provided, we need to delete it.
if (!multiplyResultProvided) {
delete copyX;
// If we allocated auxiliary memory, we need to dispose of it now.
if (allocatedAuxMemory) {
this->deallocateAuxMemory(LinearEquationSolverOperation::MultiplyRepeatedly);
}
}
template<typename ValueType>
bool LinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const {
if (!auxiliaryRepeatedMultiplyMemory) {
auxiliaryRepeatedMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixColumnCount());
return true;
}
return false;
}
template<typename ValueType>
bool LinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const {
if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) {
if (auxiliaryRepeatedMultiplyMemory) {
auxiliaryRepeatedMultiplyMemory.reset();
return true;
}
}
return false;
}
template<typename ValueType>
bool LinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) {
if (auxiliaryRepeatedMultiplyMemory) {
result = auxiliaryRepeatedMultiplyMemory->size() != this->getMatrixColumnCount();
auxiliaryRepeatedMultiplyMemory->resize(this->getMatrixColumnCount());
}
}
return result;
}
template<typename ValueType>
bool LinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const {
if (operation == LinearEquationSolverOperation::MultiplyRepeatedly) {
return static_cast<bool>(auxiliaryRepeatedMultiplyMemory);
}
return false;
}
template<typename ValueType>

72
src/solver/LinearEquationSolver.h

@ -11,6 +11,10 @@
namespace storm {
namespace solver {
enum class LinearEquationSolverOperation {
SolveEquations, MultiplyRepeatedly
};
/*!
* An interface that represents an abstract linear equation solver. In addition to solving a system of linear
* equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
@ -18,10 +22,15 @@ namespace storm {
template<class ValueType>
class LinearEquationSolver : public AbstractEquationSolver<ValueType> {
public:
LinearEquationSolver();
virtual ~LinearEquationSolver() {
// Intentionally left empty.
}
virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) = 0;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) = 0;
/*!
* Solves the equation system A*x = b. The matrix A is required to be square and have a unique solution.
* The solution of the set of linear equations will be written to the vector x. Note that the matrix A has
@ -29,22 +38,20 @@ namespace storm {
*
* @param x The solution vector that has to be computed. Its length must be equal to the number of rows of A.
* @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A.
* @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
* vector must be equal to the number of rows of A.
*/
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const = 0;
/*!
* Performs on matrix-vector multiplication x' = A*x + b.
*
* @param x The input vector with which to multiply the matrix. Its length must be equal
* to the number of columns of A.
* @param result The target vector into which to write the multiplication result. Its length must be equal
* to the number of rows of A.
* @param b If non-null, this vector is added after the multiplication. If given, its length must be equal
* to the number of rows of A.
* @param result The target vector into which to write the multiplication result. Its length must be equal
* to the number of rows of A.
*/
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const = 0;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const = 0;
/*!
* Performs repeated matrix-vector multiplication, using x[0] = x and x[i + 1] = A*x[i] + b. After
@ -55,10 +62,57 @@ namespace storm {
* to the number of columns of A.
* @param b If non-null, this vector is added after each multiplication. If given, its length must be equal
* to the number of rows of A.
* @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
* vector must be equal to the number of rows of A.
* @param n The number of times to perform the multiplication.
*/
void repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, uint_fast64_t n) const;
// Methods related to allocating/freeing auxiliary storage.
/*!
* Allocates auxiliary memory that can be used to perform the provided operation. Repeated calls to the
* corresponding function can then be run without allocating/deallocating this memory repeatedly.
* Note: Since the allocated memory is fit to the currently selected options of the solver, they must not
* be changed any more after allocating the auxiliary memory until it is deallocated again.
*
* @return True iff auxiliary memory was allocated.
*/
virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const;
/*!
* Destroys previously allocated auxiliary memory for the provided operation.
*
* @return True iff auxiliary memory was deallocated.
*/
void repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType> const* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const;
virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const;
/*!
* If the matrix dimensions changed and auxiliary memory was allocated, this function needs to be called to
* update the auxiliary memory.
*
* @return True iff the auxiliary memory was reallocated.
*/
virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const;
/*!
* Checks whether the solver has allocated auxiliary memory for the provided operation.
*
* @return True iff auxiliary memory was previously allocated (and not yet deallocated).
*/
virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const;
private:
/*!
* Retrieves the row count of the matrix associated with this solver.
*/
virtual uint64_t getMatrixRowCount() const = 0;
/*!
* Retrieves the column count of the matrix associated with this solver.
*/
virtual uint64_t getMatrixColumnCount() const = 0;
// Auxiliary memory for repeated matrix-vector multiplication.
mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRepeatedMultiplyMemory;
};
template<typename ValueType>

23
src/solver/MinMaxLinearEquationSolver.cpp

@ -28,15 +28,15 @@ namespace storm {
}
template<typename ValueType>
void MinMaxLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
void MinMaxLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
STORM_LOG_THROW(isSet(this->direction), storm::exceptions::IllegalFunctionCallException, "Optimization direction not set.");
solveEquations(convert(this->direction), x, b, multiplyResult, newX);
solveEquations(convert(this->direction), x, b);
}
template<typename ValueType>
void MinMaxLinearEquationSolver<ValueType>::multiply( std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
void MinMaxLinearEquationSolver<ValueType>::repeatedMultiply( std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const {
STORM_LOG_THROW(isSet(this->direction), storm::exceptions::IllegalFunctionCallException, "Optimization direction not set.");
return multiply(convert(this->direction), x, b, n, multiplyResult);
return repeatedMultiply(convert(this->direction), x, b, n);
}
template<typename ValueType>
@ -79,6 +79,21 @@ namespace storm {
return std::move(scheduler.get());
}
template<typename ValueType>
bool MinMaxLinearEquationSolver<ValueType>::allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
return false;
}
template<typename ValueType>
bool MinMaxLinearEquationSolver<ValueType>::deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
return false;
}
template<typename ValueType>
bool MinMaxLinearEquationSolver<ValueType>::hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
return false;
}
template<typename ValueType>
MinMaxLinearEquationSolverFactory<ValueType>::MinMaxLinearEquationSolverFactory(bool trackScheduler) : trackScheduler(trackScheduler) {
// Intentionally left empty.

45
src/solver/MinMaxLinearEquationSolver.h

@ -20,6 +20,10 @@ namespace storm {
namespace solver {
enum class MinMaxLinearEquationSolverOperation {
SolveEquations, MultiplyRepeatedly
};
/*!
* A class representing the interface that all min-max linear equation solvers shall implement.
*/
@ -40,20 +44,15 @@ namespace storm {
* @param x The solution vector x. The initial values of x represent a guess of the real values to the
* solver, but may be ignored.
* @param b The vector to add after matrix-vector multiplication.
* @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
* vector must be equal to the number of rows of A.
* @param newX If non-null, this memory is used as a scratch memory. If given, the length of this
* vector must be equal to the length of the vector x (and thus to the number of columns of A).
* @return The solution vector x of the system of linear equations as the content of the parameter x.
*/
virtual void solveEquations(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const = 0;
virtual void solveEquations(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType> const& b) const = 0;
/*!
* Behaves the same as the other variant of <code>solveEquations</code>, with the distinction that
* instead of providing the optimization direction as an argument, the internally set optimization direction
* is used. Note: this method can only be called after setting the optimization direction.
*/
void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const;
void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
/*!
* Performs (repeated) matrix-vector multiplication with the given parameters, i.e. computes
@ -68,18 +67,16 @@ namespace storm {
* i.e. after the method returns, this vector will contain the computed values.
* @param b If not null, this vector is added after each multiplication.
* @param n Specifies the number of iterations the matrix-vector multiplication is performed.
* @param multiplyResult If non-null, this memory is used as a scratch memory. If given, the length of this
* vector must be equal to the number of rows of A.
* @return The result of the repeated matrix-vector multiplication as the content of the vector x.
*/
virtual void multiply(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const = 0;
virtual void repeatedMultiply(OptimizationDirection d, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n = 1) const = 0;
/*!
* Behaves the same as the other variant of <code>multiply</code>, with the
* distinction that instead of providing the optimization direction as an argument, the internally set
* optimization direction is used. Note: this method can only be called after setting the optimization direction.
*/
virtual void multiply( std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const;
virtual void repeatedMultiply(std::vector<ValueType>& x, std::vector<ValueType>* b , uint_fast64_t n) const;
/*!
* Sets an optimization direction to use for calls to methods that do not explicitly provide one.
@ -119,6 +116,32 @@ namespace storm {
*/
std::unique_ptr<storm::storage::TotalScheduler> getScheduler();
// Methods related to allocating/freeing auxiliary storage.
/*!
* Allocates auxiliary storage that can be used to perform the provided operation. Repeated calls to the
* corresponding function can then be run without allocating/deallocating this storage repeatedly.
* Note: Since the allocated storage is fit to the currently selected options of the solver, they must not
* be changed any more after allocating the auxiliary storage until the storage is deallocated again.
*
* @return True iff auxiliary storage was allocated.
*/
virtual bool allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const;
/*!
* Destroys previously allocated auxiliary storage for the provided operation.
*
* @return True iff auxiliary storage was deallocated.
*/
virtual bool deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const;
/*!
* Checks whether the solver has allocated auxiliary storage for the provided operation.
*
* @return True iff auxiliary storage was previously allocated (and not yet deallocated).
*/
virtual bool hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const;
protected:
/// The optimization direction to use for calls to functions that do not provide it explicitly. Can also be unset.
OptimizationDirectionSetting direction;

153
src/solver/NativeLinearEquationSolver.cpp

@ -84,70 +84,62 @@ namespace storm {
}
template<typename ValueType>
NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(A), settings(settings) {
// Intentionally left empty.
NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingMemory(nullptr) {
this->setMatrix(A);
}
template<typename ValueType>
NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A))), A(*localA), settings(settings) {
// Intentionally left empty.
NativeLinearEquationSolver<ValueType>::NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings) : localA(nullptr), A(nullptr), settings(settings), auxiliarySolvingMemory(nullptr) {
this->setMatrix(std::move(A));
}
template<typename ValueType>
void NativeLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
void NativeLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType> const& A) {
localA.reset();
this->A = &A;
}
template<typename ValueType>
void NativeLinearEquationSolver<ValueType>::setMatrix(storm::storage::SparseMatrix<ValueType>&& A) {
localA = std::make_unique<storm::storage::SparseMatrix<ValueType>>(std::move(A));
this->A = localA.get();
}
template<typename ValueType>
void NativeLinearEquationSolver<ValueType>::solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
bool allocatedAuxStorage = !this->hasAuxMemory(LinearEquationSolverOperation::SolveEquations);
if (allocatedAuxStorage) {
this->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
}
if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::GaussSeidel) {
// Define the omega used for SOR.
ValueType omega = this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings<ValueType>::SolutionMethod::SOR ? this->getSettings().getOmega() : storm::utility::one<ValueType>();
// To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
bool tmpXProvided = true;
std::vector<ValueType>* tmpX = multiplyResult;
if (multiplyResult == nullptr) {
tmpX = new std::vector<ValueType>(x);
tmpXProvided = false;
} else {
*tmpX = x;
}
// Set up additional environment variables.
uint_fast64_t iterationCount = 0;
bool converged = false;
while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations()) {
A.performSuccessiveOverRelaxationStep(omega, x, b);
A->performSuccessiveOverRelaxationStep(omega, x, b);
// Now check if the process already converged within our precision.
converged = storm::utility::vector::equalModuloPrecision<ValueType>(x, *tmpX, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x));
converged = storm::utility::vector::equalModuloPrecision<ValueType>(*auxiliarySolvingMemory, x, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x));
// If we did not yet converge, we need to copy the contents of x to *tmpX.
// If we did not yet converge, we need to backup the contents of x.
if (!converged) {
*tmpX = x;
*auxiliarySolvingMemory = x;
}
// Increase iteration count so we can abort if convergence is too slow.
++iterationCount;
}
// If the vector for the temporary multiplication result was not provided, we need to delete it.
if (!tmpXProvided) {
delete tmpX;
}
} else {
// Get a Jacobi decomposition of the matrix A.
std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A->getJacobiDecomposition();
// To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
bool multiplyResultProvided = true;
std::vector<ValueType>* nextX = multiplyResult;
if (nextX == nullptr) {
nextX = new std::vector<ValueType>(x.size());
multiplyResultProvided = false;
}
std::vector<ValueType> const* copyX = nextX;
std::vector<ValueType>* currentX = &x;
// Target vector for multiplication.
std::vector<ValueType> tmpX(x.size());
std::vector<ValueType>* nextX = auxiliarySolvingMemory.get();
// Set up additional environment variables.
uint_fast64_t iterationCount = 0;
@ -155,15 +147,15 @@ namespace storm {
while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) {
// Compute D^-1 * (b - LU * x) and store result in nextX.
jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
storm::utility::vector::subtractVectors(b, tmpX, tmpX);
storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
// Swap the two pointers as a preparation for the next iteration.
std::swap(nextX, currentX);
jacobiDecomposition.first.multiplyWithVector(*currentX, *nextX);
storm::utility::vector::subtractVectors(b, *nextX, *nextX);
storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, *nextX, *nextX);
// Now check if the process already converged within our precision.
converged = storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *nextX, static_cast<ValueType>(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion());
// Swap the two pointers as a preparation for the next iteration.
std::swap(nextX, currentX);
// Increase iteration count so we can abort if convergence is too slow.
++iterationCount;
@ -171,28 +163,28 @@ namespace storm {
// If the last iteration did not write to the original x we have to swap the contents, because the
// output has to be written to the input parameter x.
if (currentX == copyX) {
if (currentX == auxiliarySolvingMemory.get()) {
std::swap(x, *currentX);
}
// If the vector for the temporary multiplication result was not provided, we need to delete it.
if (!multiplyResultProvided) {
delete copyX;
}
}
// If we allocated auxiliary memory, we need to dispose of it now.
if (allocatedAuxStorage) {
this->deallocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
}
}
template<typename ValueType>
void NativeLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b) const {
void NativeLinearEquationSolver<ValueType>::multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const {
if (&x != &result) {
A.multiplyWithVector(x, result);
A->multiplyWithVector(x, result);
if (b != nullptr) {
storm::utility::vector::addVectors(result, *b, result);
}
} else {
// If the two vectors are aliases, we need to create a temporary.
std::vector<ValueType> tmp(result.size());
A.multiplyWithVector(x, tmp);
A->multiplyWithVector(x, tmp);
if (b != nullptr) {
storm::utility::vector::addVectors(tmp, *b, result);
}
@ -209,6 +201,65 @@ namespace storm {
return settings;
}
template<typename ValueType>
bool NativeLinearEquationSolver<ValueType>::allocateAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::SolveEquations) {
if (!auxiliarySolvingMemory) {
auxiliarySolvingMemory = std::make_unique<std::vector<ValueType>>(this->getMatrixRowCount());
result = true;
}
}
result |= LinearEquationSolver<ValueType>::allocateAuxMemory(operation);
return result;
}
template<typename ValueType>
bool NativeLinearEquationSolver<ValueType>::deallocateAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::SolveEquations) {
if (auxiliarySolvingMemory) {
result = true;
auxiliarySolvingMemory.reset();
}
}
result |= LinearEquationSolver<ValueType>::deallocateAuxMemory(operation);
return result;
}
template<typename ValueType>
bool NativeLinearEquationSolver<ValueType>::reallocateAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::SolveEquations) {
if (auxiliarySolvingMemory) {
result = auxiliarySolvingMemory->size() != this->getMatrixColumnCount();
auxiliarySolvingMemory->resize(this->getMatrixRowCount());
}
}
result |= LinearEquationSolver<ValueType>::reallocateAuxMemory(operation);
return result;
}
template<typename ValueType>
bool NativeLinearEquationSolver<ValueType>::hasAuxMemory(LinearEquationSolverOperation operation) const {
bool result = false;
if (operation == LinearEquationSolverOperation::SolveEquations) {
result |= static_cast<bool>(auxiliarySolvingMemory);
}
result |= LinearEquationSolver<ValueType>::hasAuxMemory(operation);
return result;
}
template<typename ValueType>
uint64_t NativeLinearEquationSolver<ValueType>::getMatrixRowCount() const {
return this->A->getRowCount();
}
template<typename ValueType>
uint64_t NativeLinearEquationSolver<ValueType>::getMatrixColumnCount() const {
return this->A->getColumnCount();
}
template<typename ValueType>
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> NativeLinearEquationSolverFactory<ValueType>::create(storm::storage::SparseMatrix<ValueType> const& matrix) const {
return std::make_unique<storm::solver::NativeLinearEquationSolver<ValueType>>(matrix, settings);

24
src/solver/NativeLinearEquationSolver.h

@ -46,23 +46,37 @@ namespace storm {
NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, NativeLinearEquationSolverSettings<ValueType> const& settings = NativeLinearEquationSolverSettings<ValueType>());
NativeLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, NativeLinearEquationSolverSettings<ValueType> const& settings = NativeLinearEquationSolverSettings<ValueType>());
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType>& result, std::vector<ValueType> const* b = nullptr) const override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType> const& A) override;
virtual void setMatrix(storm::storage::SparseMatrix<ValueType>&& A) override;
virtual void solveEquations(std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
virtual void multiply(std::vector<ValueType>& x, std::vector<ValueType> const* b, std::vector<ValueType>& result) const override;
NativeLinearEquationSolverSettings<ValueType>& getSettings();
NativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
virtual bool allocateAuxMemory(LinearEquationSolverOperation operation) const override;
virtual bool deallocateAuxMemory(LinearEquationSolverOperation operation) const override;
virtual bool reallocateAuxMemory(LinearEquationSolverOperation operation) const override;
virtual bool hasAuxMemory(LinearEquationSolverOperation operation) const override;
private:
virtual uint64_t getMatrixRowCount() const override;
virtual uint64_t getMatrixColumnCount() const override;
// If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted
// when the solver is destructed.
std::unique_ptr<storm::storage::SparseMatrix<ValueType>> localA;
// A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
// the reference refers to localA.
storm::storage::SparseMatrix<ValueType> const& A;
// A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix
// the pointer refers to localA.
storm::storage::SparseMatrix<ValueType> const* A;
// The settings used by the solver.
NativeLinearEquationSolverSettings<ValueType> settings;
// Auxiliary memory for the equation solving methods.
mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingMemory;
};
template<typename ValueType>

172
src/solver/StandardMinMaxLinearEquationSolver.cpp

@ -84,45 +84,38 @@ namespace storm {
}
template<typename ValueType>
void StandardMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
void StandardMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
switch (this->getSettings().getSolutionMethod()) {
case StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration: solveEquationsValueIteration(dir, x, b, multiplyResult, newX); break;
case StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration: solveEquationsPolicyIteration(dir, x, b, multiplyResult, newX); break;
case StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration: solveEquationsValueIteration(dir, x, b); break;
case StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration: solveEquationsPolicyIteration(dir, x, b); break;
}
}
template<typename ValueType>
void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
// Create scratch memory if none was provided.
bool multiplyResultMemoryProvided = multiplyResult != nullptr;
if (multiplyResult == nullptr) {
multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
}
void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// Create the initial scheduler.
std::vector<storm::storage::sparse::state_type> scheduler(this->A.getRowGroupCount());
// Create a vector for storing the right-hand side of the inner equation system.
std::vector<ValueType> subB(this->A.getRowGroupCount());
// Create a vector that the inner equation solver can use as scratch memory.
std::vector<ValueType> deterministicMultiplyResult(this->A.getRowGroupCount());
// Resolve the nondeterminism according to the current scheduler.
storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(scheduler, true);
submatrix.convertToEquationSystem();
storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A.getRowGroupIndices(), b);
// Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration.
auto solver = linearEquationSolverFactory->create(std::move(submatrix));
solver->allocateAuxMemory(LinearEquationSolverOperation::SolveEquations);
Status status = Status::InProgress;
uint64_t iterations = 0;
do {
// Resolve the nondeterminism according to the current scheduler.
storm::storage::SparseMatrix<ValueType> submatrix = this->A.selectRowsFromRowGroups(scheduler, true);
submatrix.convertToEquationSystem();
storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A.getRowGroupIndices(), b);
// Solve the equation system for the 'DTMC'.
// FIXME: we need to remove the 0- and 1- states to make the solution unique.
// HOWEVER: if we start with a valid scheduler, then we will never get an illegal one, because staying
// within illegal MECs will never strictly improve the value. Is this true?
auto solver = linearEquationSolverFactory->create(std::move(submatrix));
solver->solveEquations(x, subB, &deterministicMultiplyResult);
solver->solveEquations(x, subB);
// Go through the multiplication result and see whether we can improve any of the choices.
bool schedulerImproved = false;
@ -151,6 +144,12 @@ namespace storm {
// If the scheduler did not improve, we are done.
if (!schedulerImproved) {
status = Status::Converged;
} else {
// Update the scheduler and the solver.
submatrix = this->A.selectRowsFromRowGroups(scheduler, true);
submatrix.convertToEquationSystem();
storm::utility::vector::selectVectorValues<ValueType>(subB, scheduler, this->A.getRowGroupIndices(), b);
solver->setMatrix(std::move(submatrix));
}
// Update environment variables.
@ -164,10 +163,6 @@ namespace storm {
if (this->isTrackSchedulerSet()) {
this->scheduler = std::make_unique<storm::storage::TotalScheduler>(std::move(scheduler));
}
if (!multiplyResultMemoryProvided) {
delete multiplyResult;
}
}
template<typename ValueType>
@ -186,22 +181,15 @@ namespace storm {
}
template<typename ValueType>
void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
void StandardMinMaxLinearEquationSolver<ValueType>::solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
// Create scratch memory if none was provided.
bool multiplyResultMemoryProvided = multiplyResult != nullptr;
if (multiplyResult == nullptr) {
multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
}
std::vector<ValueType>* currentX = &x;
bool xMemoryProvided = newX != nullptr;
if (newX == nullptr) {
newX = new std::vector<ValueType>(x.size());
bool allocatedAuxMemory = !this->hasAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations);
if (allocatedAuxMemory) {
this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations);
}
// Keep track of which of the vectors for x is the auxiliary copy.
std::vector<ValueType>* copyX = newX;
std::vector<ValueType>* currentX = &x;
std::vector<ValueType>* newX = auxiliarySolvingVectorMemory.get();
// Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations.
uint64_t iterations = 0;
@ -209,10 +197,10 @@ namespace storm {
Status status = Status::InProgress;
while (status == Status::InProgress) {
// Compute x' = A*x + b.
solver->multiply(*currentX, *multiplyResult, &b);
solver->multiply(*currentX, &b, *auxiliarySolvingMultiplyMemory);
// Reduce the vector x' by applying min/max for all non-deterministic choices.
storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, *newX, this->A.getRowGroupIndices());
storm::utility::vector::reduceVectorMinOrMax(dir, *auxiliarySolvingMultiplyMemory, *newX, this->A.getRowGroupIndices());
// Determine whether the method converged.
if (storm::utility::vector::equalModuloPrecision<ValueType>(*currentX, *newX, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion())) {
@ -229,10 +217,10 @@ namespace storm {
// If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result
// is currently stored in currentX, but x is the output vector.
if (currentX == copyX) {
if (currentX == auxiliarySolvingVectorMemory.get()) {
std::swap(x, *currentX);
}
// If requested, we store the scheduler for retrieval.
if (this->isTrackSchedulerSet()) {
if(iterations==0){ //may happen due to custom termination condition. Then we need to compute x'= A*x+b
@ -243,36 +231,33 @@ namespace storm {
storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices(), &choices);
this->scheduler = std::make_unique<storm::storage::TotalScheduler>(std::move(choices));
}
// Dispose of allocated scratch memory.
if (!xMemoryProvided) {
delete copyX;
}
if (!multiplyResultMemoryProvided) {
delete multiplyResult;
// If we allocated auxiliary memory, we need to dispose of it now.
if (allocatedAuxMemory) {
this->deallocateAuxMemory(MinMaxLinearEquationSolverOperation::SolveEquations);
}
}
template<typename ValueType>
void StandardMinMaxLinearEquationSolver<ValueType>::multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
// If scratch memory was not provided, we create it.
bool multiplyResultMemoryProvided = multiplyResult != nullptr;
if (!multiplyResult) {
multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
void StandardMinMaxLinearEquationSolver<ValueType>::repeatedMultiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const {
bool allocatedAuxMemory = !this->hasAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly);
if (allocatedAuxMemory) {
this->allocateAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly);
}
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory->create(A);
for (uint64_t i = 0; i < n; ++i) {
solver->multiply(x, *multiplyResult, b);
solver->multiply(x, b, *auxiliaryRepeatedMultiplyMemory);
// Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
// element of the min/max operator stack.
storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices());
storm::utility::vector::reduceVectorMinOrMax(dir, *auxiliaryRepeatedMultiplyMemory, x, this->A.getRowGroupIndices());
}
if (!multiplyResultMemoryProvided) {
delete multiplyResult;
// If we allocated auxiliary memory, we need to dispose of it now.
if (allocatedAuxMemory) {
this->deallocateAuxMemory(MinMaxLinearEquationSolverOperation::MultiplyRepeatedly);
}
}
@ -309,6 +294,73 @@ namespace storm {
return settings;
}
template<typename ValueType>
bool StandardMinMaxLinearEquationSolver<ValueType>::allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
bool result = false;
if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) {
if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) {
if (!auxiliarySolvingMultiplyMemory) {
result = true;
auxiliarySolvingMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowCount());
}
if (!auxiliarySolvingVectorMemory) {
result = true;
auxiliarySolvingVectorMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowGroupCount());
}
} else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) {
// Nothing to do in this case.
} else {
STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method.");
}
} else {
if (!auxiliaryRepeatedMultiplyMemory) {
result = true;
auxiliaryRepeatedMultiplyMemory = std::make_unique<std::vector<ValueType>>(this->A.getRowCount());
}
}
return result;
}
template<typename ValueType>
bool StandardMinMaxLinearEquationSolver<ValueType>::deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
bool result = false;
if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) {
if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) {
if (auxiliarySolvingMultiplyMemory) {
result = true;
auxiliarySolvingMultiplyMemory.reset();
}
if (auxiliarySolvingVectorMemory) {
result = true;
auxiliarySolvingVectorMemory.reset();
}
} else if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::PolicyIteration) {
// Nothing to do in this case.
} else {
STORM_LOG_ASSERT(false, "Cannot allocate aux storage for this method.");
}
} else {
if (auxiliaryRepeatedMultiplyMemory) {
result = true;
auxiliaryRepeatedMultiplyMemory.reset();
}
}
return result;
}
template<typename ValueType>
bool StandardMinMaxLinearEquationSolver<ValueType>::hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const {
if (operation == MinMaxLinearEquationSolverOperation::SolveEquations) {
if (this->getSettings().getSolutionMethod() == StandardMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) {
return auxiliarySolvingMultiplyMemory && auxiliarySolvingVectorMemory;
} else {
return false;
}
} else {
return static_cast<bool>(auxiliaryRepeatedMultiplyMemory);
}
}
template<typename ValueType>
StandardMinMaxLinearEquationSolverFactory<ValueType>::StandardMinMaxLinearEquationSolverFactory(bool trackScheduler) : MinMaxLinearEquationSolverFactory<ValueType>(trackScheduler), linearEquationSolverFactory(nullptr) {
// Intentionally left empty.

19
src/solver/StandardMinMaxLinearEquationSolver.h

@ -38,15 +38,19 @@ namespace storm {
StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings<ValueType> const& settings = StandardMinMaxLinearEquationSolverSettings<ValueType>());
StandardMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType>&& A, std::unique_ptr<LinearEquationSolverFactory<ValueType>>&& linearEquationSolverFactory, StandardMinMaxLinearEquationSolverSettings<ValueType> const& settings = StandardMinMaxLinearEquationSolverSettings<ValueType>());
virtual void solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
virtual void multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b = nullptr, uint_fast64_t n = 1, std::vector<ValueType>* multiplyResult = nullptr) const override;
virtual void solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
virtual void repeatedMultiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const override;
StandardMinMaxLinearEquationSolverSettings<ValueType> const& getSettings() const;
StandardMinMaxLinearEquationSolverSettings<ValueType>& getSettings();
virtual bool allocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const override;
virtual bool deallocateAuxMemory(MinMaxLinearEquationSolverOperation operation) const override;
virtual bool hasAuxMemory(MinMaxLinearEquationSolverOperation operation) const override;
private:
void solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const;
void solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const;
void solveEquationsPolicyIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
void solveEquationsValueIteration(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const;
@ -70,6 +74,13 @@ namespace storm {
// A reference to the original sparse matrix given to this solver. If the solver takes posession of the matrix
// the reference refers to localA.
storm::storage::SparseMatrix<ValueType> const& A;
// Auxiliary memory for equation solving.
mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingMultiplyMemory;
mutable std::unique_ptr<std::vector<ValueType>> auxiliarySolvingVectorMemory;
// Auxiliary memory for repeated matrix-vector multiplication.
mutable std::unique_ptr<std::vector<ValueType>> auxiliaryRepeatedMultiplyMemory;
};
template<typename ValueType>

18
src/solver/TopologicalMinMaxLinearEquationSolver.cpp

@ -33,7 +33,7 @@ namespace storm {
}
template<typename ValueType>
void TopologicalMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult, std::vector<ValueType>* newX) const {
void TopologicalMinMaxLinearEquationSolver<ValueType>::solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
#ifdef GPU_USE_FLOAT
#define __FORCE_FLOAT_CALCULATION true
@ -49,7 +49,7 @@ namespace storm {
std::vector<float> new_x = storm::utility::vector::toValueType<float>(x);
std::vector<float> const new_b = storm::utility::vector::toValueType<float>(b);
newSolver.solveEquations(dir, new_x, new_b, nullptr, nullptr);
newSolver.solveEquations(dir, new_x, new_b);
for (size_t i = 0, size = new_x.size(); i < size; ++i) {
x.at(i) = new_x.at(i);
@ -422,14 +422,8 @@ namespace storm {
}
template<typename ValueType>
void TopologicalMinMaxLinearEquationSolver<ValueType>::multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const {
// If scratch memory was not provided, we need to create it.
bool multiplyResultMemoryProvided = true;
if (multiplyResult == nullptr) {
multiplyResult = new std::vector<ValueType>(this->A.getRowCount());
multiplyResultMemoryProvided = false;
}
void TopologicalMinMaxLinearEquationSolver<ValueType>::repeatedMultiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const {
std::unique_ptr<std::vector<ValueType>> multiplyResult = std::make_unique<std::vector<ValueType>>(this->A.getRowCount());
// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
for (uint_fast64_t i = 0; i < n; ++i) {
@ -445,10 +439,6 @@ namespace storm {
storm::utility::vector::reduceVectorMinOrMax(dir, *multiplyResult, x, this->A.getRowGroupIndices());
}
if (!multiplyResultMemoryProvided) {
delete multiplyResult;
}
}
template<typename ValueType>

4
src/solver/TopologicalMinMaxLinearEquationSolver.h

@ -32,9 +32,9 @@ namespace storm {
*/
TopologicalMinMaxLinearEquationSolver(storm::storage::SparseMatrix<ValueType> const& A, double precision = 1e-6, uint_fast64_t maximalNumberOfIterations = 20000, bool relative = true);
virtual void solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult = nullptr, std::vector<ValueType>* newX = nullptr) const override;
virtual void solveEquations(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const override;
virtual void multiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n, std::vector<ValueType>* multiplyResult) const override;
virtual void repeatedMultiply(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType>* b, uint_fast64_t n) const override;
private:
storm::storage::SparseMatrix<ValueType> const& A;

10
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.");

11
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.

64
src/storage/SparseMatrix.cpp

@ -658,6 +658,70 @@ namespace storm {
return bv;
}
template<typename ValueType>
void SparseMatrix<ValueType>::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<MatrixEntry<index_type, value_type>> 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<typename ValueType>
ValueType SparseMatrix<ValueType>::getConstrainedRowSum(index_type row, storm::storage::BitVector const& constraint) const {
ValueType result = storm::utility::zero<ValueType>();

8
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.
*

2
src/utility/ModelInstantiator.h

@ -82,7 +82,7 @@ namespace storm {
auto markovianStatesCopy = parametricModel.getMarkovianStates();
auto choiceLabelingCopy = parametricModel.getOptionalChoiceLabeling();
std::vector<ConstantType> exitRates(parametricModel.getExitRates().size(), storm::utility::one<ConstantType>());
this->instantiatedModel = std::make_shared<ConstantSparseModelType>(buildDummyMatrix(parametricModel.getTransitionMatrix()), std::move(stateLabelingCopy), std::move(markovianStatesCopy), std::move(exitRates), buildDummyRewardModels(parametricModel.getRewardModels()), std::move(choiceLabelingCopy));
this->instantiatedModel = std::make_shared<ConstantSparseModelType>(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());
}

26
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<storm::models::sparse::Model<double>> model = storm::builder::ExplicitModelBuilder<double>(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<storm::models::sparse::MarkovAutomaton<double>>()->getMarkovianStates().getNumberOfSetBits());
program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/hybrid_states.ma");
model = storm::builder::ExplicitModelBuilder<double>(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<storm::models::sparse::MarkovAutomaton<double>>()->getMarkovianStates().getNumberOfSetBits());
program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/stream2.ma");
model = storm::builder::ExplicitModelBuilder<double>(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<storm::models::sparse::MarkovAutomaton<double>>()->getMarkovianStates().getNumberOfSetBits());
}
TEST(ExplicitPrismModelBuilderTest, FailComposition) {
storm::prism::Program program = storm::parser::PrismParser::parse(STORM_CPP_TESTS_BASE_PATH "/functional/builder/system_composition.nm");

19
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

15
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

50
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<N -> 1 : (s'=1);
[start] s=0 & k<n -> 1 : (s'=2) & (k'=k+1);
<> s=1 -> inRate : (n'=n+1) & (s'=0);
<> s=2 & n<N & k<n -> inRate : (n'=n+1) + processingRate : (k'=k+1);
<> s=2 & n<N & k=n -> inRate : (n'=n+1) + processingRate : (s'=0);
<> s=2 & n=N & k<n -> 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

10
test/functional/solver/GmmxxMinMaxLinearEquationSolverTest.cpp

@ -78,23 +78,23 @@ TEST(GmmxxMinMaxLinearEquationSolver, MatrixVectorMultiplication) {
auto factory = storm::solver::GmmxxMinMaxLinearEquationSolverFactory<double>();
auto solver = factory.create(A);
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 1));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 1));
ASSERT_LT(std::abs(x[0] - 0.099), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
x = {0, 1, 0};
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 2));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 2));
ASSERT_LT(std::abs(x[0] - 0.1881), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
x = {0, 1, 0};
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 20));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 20));
ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
x = {0, 1, 0};
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Maximize, x, nullptr, 1));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Maximize, x, nullptr, 1));
ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
x = {0, 1, 0};
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Maximize, x, nullptr, 20));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Maximize, x, nullptr, 20));
ASSERT_LT(std::abs(x[0] - 0.9238082658), storm::settings::getModule<storm::settings::modules::GmmxxEquationSolverSettings>().getPrecision());
}

10
test/functional/solver/NativeMinMaxLinearEquationSolverTest.cpp

@ -49,23 +49,23 @@ TEST(NativeMinMaxLinearEquationSolver, MatrixVectorMultiplication) {
auto factory = storm::solver::NativeMinMaxLinearEquationSolverFactory<double>();
auto solver = factory.create(A);
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 1));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 1));
ASSERT_LT(std::abs(x[0] - 0.099), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
x = {0, 1, 0};
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 2));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 2));
ASSERT_LT(std::abs(x[0] - 0.1881), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
x = {0, 1, 0};
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Minimize, x, nullptr, 20));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Minimize, x, nullptr, 20));
ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
x = {0, 1, 0};
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Maximize, x, nullptr, 1));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Maximize, x, nullptr, 1));
ASSERT_LT(std::abs(x[0] - 0.5), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
x = {0, 1, 0};
ASSERT_NO_THROW(solver->multiply(storm::OptimizationDirection::Maximize, x, nullptr, 20));
ASSERT_NO_THROW(solver->repeatedMultiply(storm::OptimizationDirection::Maximize, x, nullptr, 20));
ASSERT_LT(std::abs(x[0] - 0.9238082658), storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>().getPrecision());
}

Loading…
Cancel
Save