From 3b4dca1a03419579583e15ee2ff398bc80a46f3f Mon Sep 17 00:00:00 2001
From: dehnert <dehnert@cs.rwth-aachen.de>
Date: Fri, 3 Apr 2015 16:12:11 +0200
Subject: [PATCH] Improved Jacobi method a bit.
Former-commit-id: f4affeebf67891dfca5930abf2dcc915a396e4e9
---
.../csl/SparseCtmcCslModelChecker.cpp | 10 +--
.../SparseMarkovAutomatonCslModelChecker.cpp | 10 +--
.../prctl/HybridDtmcPrctlModelChecker.cpp | 10 ++-
.../prctl/SparseDtmcPrctlModelChecker.cpp | 4 +-
.../prctl/SparseMdpPrctlModelChecker.cpp | 4 +-
.../SparseDtmcEliminationModelChecker.cpp | 2 +-
.../SymbolicQuantitativeCheckResult.cpp | 18 ++--
src/models/sparse/Model.cpp | 2 +-
src/solver/GmmxxLinearEquationSolver.cpp | 6 +-
src/solver/NativeLinearEquationSolver.cpp | 11 +--
.../NativeMinMaxLinearEquationSolver.cpp | 4 +-
.../TopologicalMinMaxLinearEquationSolver.cpp | 2 +-
src/storage/SparseMatrix.cpp | 32 +++-----
src/storage/SparseMatrix.h | 4 +-
src/storage/dd/CuddAdd.cpp | 82 +++++++++++++++++--
src/storage/dd/CuddAdd.h | 37 +++++++++
src/storage/dd/CuddDd.cpp | 8 +-
src/storage/dd/CuddDd.h | 11 ++-
src/storage/dd/CuddDdManager.cpp | 4 +
src/storage/dd/CuddDdManager.h | 7 ++
src/storage/dd/CuddOdd.cpp | 4 +
src/storage/dd/CuddOdd.h | 7 ++
src/utility/cli.h | 2 +-
src/utility/vector.h | 63 ++++++++------
test/functional/storage/CuddDdTest.cpp | 3 +
test/functional/storage/SparseMatrixTest.cpp | 11 +--
26 files changed, 255 insertions(+), 103 deletions(-)
diff --git a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
index 412c22625..4f24cde22 100644
--- a/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseCtmcCslModelChecker.cpp
@@ -300,14 +300,14 @@ namespace storm {
storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]);
std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&foxGlynnResult] (ValueType const& a, ValueType const& b) { return a + std::get<3>(foxGlynnResult)[0] * b; };
if (addVector != nullptr) {
- storm::utility::vector::applyPointwiseInPlace(result, *addVector, addAndScale);
+ storm::utility::vector::applyPointwise(result, *addVector, result, addAndScale);
}
++startingIteration;
} else {
if (computeCumulativeReward) {
result = std::vector<ValueType>(values.size());
std::function<ValueType (ValueType const&)> scaleWithUniformizationRate = [&uniformizationRate] (ValueType const& a) -> ValueType { return a / uniformizationRate; };
- storm::utility::vector::applyPointwiseInPlace(result, scaleWithUniformizationRate);
+ storm::utility::vector::applyPointwise(result, result, scaleWithUniformizationRate);
} else {
result = std::vector<ValueType>(values.size());
}
@@ -325,7 +325,7 @@ namespace storm {
// 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->performMatrixVectorMultiplication(values, nullptr, 1, &multiplicationResult);
- storm::utility::vector::applyPointwiseInPlace(result, values, addAndScale);
+ storm::utility::vector::applyPointwise(result, values, result, addAndScale);
}
}
@@ -337,7 +337,7 @@ namespace storm {
solver->performMatrixVectorMultiplication(values, addVector, 1, &multiplicationResult);
weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
- storm::utility::vector::applyPointwiseInPlace(result, values, addAndScale);
+ storm::utility::vector::applyPointwise(result, values, result, addAndScale);
}
return result;
@@ -409,7 +409,7 @@ namespace storm {
if (this->getModel().hasTransitionRewards()) {
totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
if (this->getModel().hasStateRewards()) {
- storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector());
+ storm::utility::vector::addVectors(totalRewardVector, this->getModel().getStateRewardVector(), totalRewardVector);
}
} else {
totalRewardVector = std::vector<ValueType>(this->getModel().getStateRewardVector());
diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
index 04b36d823..abec6a929 100644
--- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
+++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.cpp
@@ -103,22 +103,22 @@ namespace storm {
for (uint_fast64_t currentStep = 0; currentStep < numberOfSteps; ++currentStep) {
// Start by (re-)computing bProbabilistic = bProbabilisticFixed + aProbabilisticToMarkovian * vMarkovian.
aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
- storm::utility::vector::addVectorsInPlace(bProbabilistic, bProbabilisticFixed);
+ storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
// Now perform the inner value iteration for probabilistic states.
solver->solveEquationSystem(min, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
// (Re-)compute bMarkovian = bMarkovianFixed + aMarkovianToProbabilistic * vProbabilistic.
aMarkovianToProbabilistic.multiplyWithVector(probabilisticNonGoalValues, bMarkovian);
- storm::utility::vector::addVectorsInPlace(bMarkovian, bMarkovianFixed);
+ storm::utility::vector::addVectors(bMarkovian, bMarkovianFixed, bMarkovian);
aMarkovian.multiplyWithVector(markovianNonGoalValues, markovianNonGoalValuesSwap);
std::swap(markovianNonGoalValues, markovianNonGoalValuesSwap);
- storm::utility::vector::addVectorsInPlace(markovianNonGoalValues, bMarkovian);
+ storm::utility::vector::addVectors(markovianNonGoalValues, bMarkovian, markovianNonGoalValues);
}
// After the loop, perform one more step of the value iteration for PS states.
aProbabilisticToMarkovian.multiplyWithVector(markovianNonGoalValues, bProbabilistic);
- storm::utility::vector::addVectorsInPlace(bProbabilistic, bProbabilisticFixed);
+ storm::utility::vector::addVectors(bProbabilistic, bProbabilisticFixed, bProbabilistic);
solver->solveEquationSystem(min, probabilisticNonGoalValues, bProbabilistic, &multiplicationResultScratchMemory, &aProbabilisticScratchMemory);
}
@@ -217,7 +217,7 @@ namespace storm {
if (model.hasTransitionRewards()) {
totalRewardVector = model.getTransitionMatrix().getPointwiseProductRowSumVector(model.getTransitionRewardMatrix());
if (model.hasStateRewards()) {
- storm::utility::vector::addVectorsInPlace(totalRewardVector, model.getStateRewardVector());
+ storm::utility::vector::addVectors(totalRewardVector, model.getStateRewardVector(), totalRewardVector);
}
} else {
totalRewardVector = std::vector<ValueType>(model.getStateRewardVector());
diff --git a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
index d71d206a4..b927031e0 100644
--- a/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/HybridDtmcPrctlModelChecker.cpp
@@ -69,18 +69,26 @@ namespace storm {
submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
+ submatrix.exportToDot("submatrix.dot");
+
// Create the solution vector.
std::vector<ValueType> x(maybeStates.getNonZeroCount(), ValueType(0.5));
// Translate the symbolic matrix/vector to their explicit representations and solve the equation system.
+ STORM_LOG_DEBUG("Converting the symbolic matrix to a sparse matrix.");
storm::storage::SparseMatrix<ValueType> explicitSubmatrix = submatrix.toMatrix(odd, odd);
+
+ STORM_LOG_DEBUG("Converting the symbolic vector to an explicit vector.");
std::vector<ValueType> b = subvector.template toVector<ValueType>(odd);
+ STORM_LOG_DEBUG("Solving explicit linear equation system.");
std::unique_ptr<storm::solver::LinearEquationSolver<ValueType>> solver = linearEquationSolverFactory.create(explicitSubmatrix);
solver->solveEquationSystem(x, b);
// Now that we have the explicit solution of the system, we need to transform it to a symbolic representation.
- storm::dd::Add<DdType> numericResult; // = storm::dd::Add<DdType>(x, odd);
+ STORM_LOG_DEBUG("Converting the explicit result to a symbolic form.");
+ storm::dd::Add<DdType> numericResult(model.getManager().asSharedPointer(), x, odd, model.getRowVariables());
+
return statesWithProbability01.second.toAdd() + numericResult;
} else {
return statesWithProbability01.second.toAdd();
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
index 1b543b5fc..46ba68183 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp
@@ -171,7 +171,7 @@ namespace storm {
if (this->getModel().hasTransitionRewards()) {
totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
if (this->getModel().hasStateRewards()) {
- storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector());
+ storm::utility::vector::addVectors(totalRewardVector, this->getModel().getStateRewardVector(), totalRewardVector);
}
} else {
totalRewardVector = std::vector<ValueType>(this->getModel().getStateRewardVector());
@@ -272,7 +272,7 @@ namespace storm {
// first.
std::vector<ValueType> subStateRewards(b.size());
storm::utility::vector::selectVectorValues(subStateRewards, maybeStates, stateRewardVector.get());
- storm::utility::vector::addVectorsInPlace(b, subStateRewards);
+ storm::utility::vector::addVectors(b, subStateRewards, b);
}
} else {
// If only a state-based reward model is available, we take this vector as the
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
index 516090698..664a15fdd 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.cpp
@@ -183,7 +183,7 @@ namespace storm {
if (this->getModel().hasTransitionRewards()) {
totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
if (this->getModel().hasStateRewards()) {
- storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector());
+ storm::utility::vector::addVectors(totalRewardVector, this->getModel().getStateRewardVector(), totalRewardVector);
}
} else {
totalRewardVector = std::vector<ValueType>(this->getModel().getStateRewardVector());
@@ -285,7 +285,7 @@ namespace storm {
// first.
std::vector<ValueType> subStateRewards(b.size());
storm::utility::vector::selectVectorValuesRepeatedly(subStateRewards, maybeStates, transitionMatrix.getRowGroupIndices(), stateRewardVector.get());
- storm::utility::vector::addVectorsInPlace(b, subStateRewards);
+ storm::utility::vector::addVectors(b, subStateRewards, b);
}
} else {
// If only a state-based reward model is available, we take this vector as the
diff --git a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
index 72dfe070e..e4087ad52 100644
--- a/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
+++ b/src/modelchecker/reachability/SparseDtmcEliminationModelChecker.cpp
@@ -176,7 +176,7 @@ namespace storm {
// first.
std::vector<ValueType> subStateRewards(stateRewards.size());
storm::utility::vector::selectVectorValues(subStateRewards, maybeStates, model.getStateRewardVector());
- storm::utility::vector::addVectorsInPlace(stateRewards, subStateRewards);
+ storm::utility::vector::addVectors(stateRewards, subStateRewards, stateRewards);
}
} else {
// If only a state-based reward model is available, we take this vector as the
diff --git a/src/modelchecker/results/SymbolicQuantitativeCheckResult.cpp b/src/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
index c265efd5c..fc18df5c6 100644
--- a/src/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
+++ b/src/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
@@ -49,14 +49,18 @@ namespace storm {
template<storm::dd::DdType Type>
std::ostream& SymbolicQuantitativeCheckResult<Type>::writeToStream(std::ostream& out) const {
out << "[";
- bool first = true;
- for (auto valuationValuePair : this->values) {
- if (!first) {
- out << ", ";
- } else {
- first = false;
+ if (this->values.isZero()) {
+ out << "0";
+ } else {
+ bool first = true;
+ for (auto valuationValuePair : this->values) {
+ if (!first) {
+ out << ", ";
+ } else {
+ first = false;
+ }
+ out << valuationValuePair.second;
}
- out << valuationValuePair.second;
}
out << "]";
return out;
diff --git a/src/models/sparse/Model.cpp b/src/models/sparse/Model.cpp
index 6541d517f..e1230d051 100644
--- a/src/models/sparse/Model.cpp
+++ b/src/models/sparse/Model.cpp
@@ -142,7 +142,7 @@ namespace storm {
void Model<ValueType>::convertTransitionRewardsToStateRewards() {
STORM_LOG_THROW(this->hasTransitionRewards(), storm::exceptions::InvalidOperationException, "Cannot reduce non-existant transition rewards to state rewards.");
if (this->hasStateRewards()) {
- storm::utility::vector::addVectorsInPlace(stateRewardVector.get(), transitionMatrix.getPointwiseProductRowSumVector(transitionRewardMatrix.get()));
+ storm::utility::vector::addVectors(stateRewardVector.get(), transitionMatrix.getPointwiseProductRowSumVector(transitionRewardMatrix.get()), stateRewardVector.get());
} else {
this->stateRewardVector = transitionMatrix.getPointwiseProductRowSumVector(transitionRewardMatrix.get());
}
diff --git a/src/solver/GmmxxLinearEquationSolver.cpp b/src/solver/GmmxxLinearEquationSolver.cpp
index 0bc5a3c46..789ee6e55 100644
--- a/src/solver/GmmxxLinearEquationSolver.cpp
+++ b/src/solver/GmmxxLinearEquationSolver.cpp
@@ -148,10 +148,8 @@ 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 {
// Get a Jacobi decomposition of the matrix A.
- std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+ std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
- // Convert the (inverted) diagonal matrix to gmm++'s format.
- std::unique_ptr<gmm::csr_matrix<ValueType>> gmmDinv = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<ValueType>(std::move(jacobiDecomposition.second));
// 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));
@@ -176,7 +174,7 @@ namespace storm {
// 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);
- gmm::mult(*gmmDinv, tmpX, *nextX);
+ storm::utility::vector::multiplyVectorsPointwise(jacobiDecomposition.second, tmpX, *nextX);
// Now check if the process already converged within our precision.
converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, precision, relative);
diff --git a/src/solver/NativeLinearEquationSolver.cpp b/src/solver/NativeLinearEquationSolver.cpp
index 47c72dae8..9206ae301 100644
--- a/src/solver/NativeLinearEquationSolver.cpp
+++ b/src/solver/NativeLinearEquationSolver.cpp
@@ -34,7 +34,9 @@ namespace storm {
template<typename ValueType>
void NativeLinearEquationSolver<ValueType>::solveEquationSystem(std::vector<ValueType>& x, std::vector<ValueType> const& b, std::vector<ValueType>* multiplyResult) const {
// Get a Jacobi decomposition of the matrix A.
- std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+ std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> jacobiDecomposition = A.getJacobiDecomposition();
+
+ std::cout << "LU has " << jacobiDecomposition.first.getNonzeroEntryCount() << " nonzeros." << std::endl;
// To avoid copying the contents of the vector in the loop, we create a temporary x to swap with.
bool multiplyResultProvided = true;
@@ -56,9 +58,8 @@ namespace storm {
while (!converged && iterationCount < maximalNumberOfIterations) {
// Compute D^-1 * (b - LU * x) and store result in nextX.
jacobiDecomposition.first.multiplyWithVector(*currentX, tmpX);
- storm::utility::vector::scaleVectorInPlace(tmpX, -storm::utility::one<ValueType>());
- storm::utility::vector::addVectorsInPlace(tmpX, b);
- jacobiDecomposition.second.multiplyWithVector(tmpX, *nextX);
+ 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);
@@ -103,7 +104,7 @@ namespace storm {
// If requested, add an offset to the current result vector.
if (b != nullptr) {
- storm::utility::vector::addVectorsInPlace(*currentX, *b);
+ storm::utility::vector::addVectors(*currentX, *b, *currentX);
}
}
diff --git a/src/solver/NativeMinMaxLinearEquationSolver.cpp b/src/solver/NativeMinMaxLinearEquationSolver.cpp
index 48ef291c7..51da3e692 100644
--- a/src/solver/NativeMinMaxLinearEquationSolver.cpp
+++ b/src/solver/NativeMinMaxLinearEquationSolver.cpp
@@ -50,7 +50,7 @@ namespace storm {
while (!converged && iterations < maximalNumberOfIterations) {
// Compute x' = A*x + b.
A.multiplyWithVector(*currentX, *multiplyResult);
- storm::utility::vector::addVectorsInPlace(*multiplyResult, b);
+ storm::utility::vector::addVectors(*multiplyResult, b, *multiplyResult);
// 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.
@@ -106,7 +106,7 @@ namespace storm {
// Add b if it is non-null.
if (b != nullptr) {
- storm::utility::vector::addVectorsInPlace(*multiplyResult, *b);
+ storm::utility::vector::addVectors(*multiplyResult, *b, *multiplyResult);
}
// Reduce the vector x' by applying min/max for all non-deterministic choices as given by the topmost
diff --git a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
index 66de10ed9..8ea3478b0 100644
--- a/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
+++ b/src/solver/TopologicalMinMaxLinearEquationSolver.cpp
@@ -254,7 +254,7 @@ namespace storm {
while (!converged && localIterations < this->maximalNumberOfIterations) {
// Compute x' = A*x + b.
sccSubmatrix.multiplyWithVector(*currentX, sccMultiplyResult);
- storm::utility::vector::addVectorsInPlace<ValueType>(sccMultiplyResult, sccSubB);
+ storm::utility::vector::addVectors<ValueType>(sccMultiplyResult, sccSubB, sccMultiplyResult);
//A.multiplyWithVector(scc, nondeterministicChoiceIndices, *currentX, multiplyResult);
//storm::utility::addVectors(scc, nondeterministicChoiceIndices, multiplyResult, b);
diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp
index 335e79910..9a8affbd1 100644
--- a/src/storage/SparseMatrix.cpp
+++ b/src/storage/SparseMatrix.cpp
@@ -11,6 +11,7 @@
#include "src/adapters/CarlAdapter.h"
#include "src/exceptions/InvalidStateException.h"
+#include "src/exceptions/NotImplementedException.h"
#include "src/utility/macros.h"
#include "log4cplus/logger.h"
@@ -701,40 +702,31 @@ namespace storm {
}
template<typename ValueType>
- typename std::pair<storm::storage::SparseMatrix<ValueType>, storm::storage::SparseMatrix<ValueType>> SparseMatrix<ValueType>::getJacobiDecomposition() const {
- if (rowCount != columnCount) {
- throw storm::exceptions::InvalidArgumentException() << "Illegal call to SparseMatrix::invertDiagonal: matrix is non-square.";
- }
- storm::storage::SparseMatrix<ValueType> resultLU(*this);
- resultLU.deleteDiagonalEntries();
+ typename std::pair<storm::storage::SparseMatrix<ValueType>, std::vector<ValueType>> SparseMatrix<ValueType>::getJacobiDecomposition() const {
+ STORM_LOG_THROW(this->getRowCount() == this->getColumnCount(), storm::exceptions::InvalidArgumentException, "Canno compute Jacobi decomposition of non-square matrix.");
- SparseMatrixBuilder<ValueType> dInvBuilder(rowCount, columnCount, rowCount);
+ // Prepare the resulting data structures.
+ SparseMatrixBuilder<ValueType> luBuilder(this->getRowCount(), this->getColumnCount());
+ std::vector<ValueType> invertedDiagonal(rowCount);
// Copy entries to the appropriate matrices.
for (index_type rowNumber = 0; rowNumber < rowCount; ++rowNumber) {
-
- // Because the matrix may have several entries on the diagonal, we need to sum them before we are able
- // to invert the entry.
- ValueType diagonalValue = storm::utility::zero<ValueType>();
for (const_iterator it = this->begin(rowNumber), ite = this->end(rowNumber); it != ite; ++it) {
if (it->getColumn() == rowNumber) {
- diagonalValue += it->getValue();
- } else if (it->getColumn() > rowNumber) {
- break;
+ invertedDiagonal[rowNumber] = storm::utility::one<ValueType>() / it->getValue();
+ } else {
+ luBuilder.addNextValue(rowNumber, it->getColumn(), it->getValue());
}
}
- dInvBuilder.addNextValue(rowNumber, rowNumber, storm::utility::one<ValueType>() / diagonalValue);
}
- return std::make_pair(std::move(resultLU), dInvBuilder.build());
+ return std::make_pair(luBuilder.build(), std::move(invertedDiagonal));
}
#ifdef STORM_HAVE_CARL
template<>
- typename std::pair<storm::storage::SparseMatrix<RationalFunction>, storm::storage::SparseMatrix<RationalFunction>> SparseMatrix<RationalFunction>::getJacobiDecomposition() const {
- // NOT SUPPORTED
- // TODO do whatever storm does in such cases.
- assert(false);
+ typename std::pair<storm::storage::SparseMatrix<RationalFunction>, std::vector<RationalFunction>> SparseMatrix<RationalFunction>::getJacobiDecomposition() const {
+ STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported.");
}
#endif
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 1bac30a6c..8eded966b 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -605,9 +605,9 @@ namespace storm {
/*!
* Calculates the Jacobi decomposition of this sparse matrix. For this operation, the matrix must be square.
*
- * @return A pair (L+U, D^-1) containing the matrix L+U and the inverted diagonal matrix D^-1.
+ * @return A pair (L+U, D^-1) containing the matrix L+U and the inverted diagonal D^-1 (as a vector).
*/
- std::pair<storm::storage::SparseMatrix<value_type>, storm::storage::SparseMatrix<value_type>> getJacobiDecomposition() const;
+ std::pair<storm::storage::SparseMatrix<value_type>, std::vector<value_type>> getJacobiDecomposition() const;
/*!
* Performs a pointwise matrix multiplication of the matrix with the given matrix and returns a vector
diff --git a/src/storage/dd/CuddAdd.cpp b/src/storage/dd/CuddAdd.cpp
index 72787075b..71ec05758 100644
--- a/src/storage/dd/CuddAdd.cpp
+++ b/src/storage/dd/CuddAdd.cpp
@@ -18,6 +18,10 @@ namespace storm {
// Intentionally left empty.
}
+ Add<DdType::CUDD>::Add(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, std::vector<double> const& values, Odd<DdType::CUDD> const& odd, std::set<storm::expressions::Variable> const& metaVariables) : Dd<DdType::CUDD>(ddManager, metaVariables) {
+ cuddAdd = fromVector(ddManager, values, odd, metaVariables);
+ }
+
Bdd<DdType::CUDD> Add<DdType::CUDD>::toBdd() const {
return Bdd<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().BddPattern(), this->getContainedMetaVariables());
}
@@ -155,19 +159,19 @@ namespace storm {
std::set_union(this->getContainedMetaVariables().begin(), this->getContainedMetaVariables().end(), other.getContainedMetaVariables().begin(), other.getContainedMetaVariables().end(), std::inserter(metaVariables, metaVariables.begin()));
return Add<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().Pow(other.getCuddAdd()), metaVariables);
}
-
+
Add<DdType::CUDD> Add<DdType::CUDD>::mod(Add<DdType::CUDD> const& other) const {
std::set<storm::expressions::Variable> metaVariables;
std::set_union(this->getContainedMetaVariables().begin(), this->getContainedMetaVariables().end(), other.getContainedMetaVariables().begin(), other.getContainedMetaVariables().end(), std::inserter(metaVariables, metaVariables.begin()));
return Add<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().Mod(other.getCuddAdd()), metaVariables);
}
-
+
Add<DdType::CUDD> Add<DdType::CUDD>::logxy(Add<DdType::CUDD> const& other) const {
std::set<storm::expressions::Variable> metaVariables;
std::set_union(this->getContainedMetaVariables().begin(), this->getContainedMetaVariables().end(), other.getContainedMetaVariables().begin(), other.getContainedMetaVariables().end(), std::inserter(metaVariables, metaVariables.begin()));
return Add<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().LogXY(other.getCuddAdd()), metaVariables);
}
-
+
Add<DdType::CUDD> Add<DdType::CUDD>::floor() const {
return Add<DdType::CUDD>(this->getDdManager(), this->getCuddAdd().Floor(), this->getContainedMetaVariables());
}
@@ -464,17 +468,20 @@ namespace storm {
ddRowVariableIndices.push_back(ddVariable.getIndex());
}
}
+ std::sort(ddRowVariableIndices.begin(), ddRowVariableIndices.end());
+
for (auto const& variable : columnMetaVariables) {
DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
for (auto const& ddVariable : metaVariable.getDdVariables()) {
ddColumnVariableIndices.push_back(ddVariable.getIndex());
}
}
+ std::sort(ddColumnVariableIndices.begin(), ddColumnVariableIndices.end());
// Prepare the vectors that represent the matrix.
std::vector<uint_fast64_t> rowIndications(rowOdd.getTotalOffset() + 1);
std::vector<storm::storage::MatrixEntry<uint_fast64_t, double>> columnsAndValues(this->getNonZeroCount());
-
+
// Create a trivial row grouping.
std::vector<uint_fast64_t> trivialRowGroupIndices(rowIndications.size());
uint_fast64_t i = 0;
@@ -526,6 +533,7 @@ namespace storm {
}
rowAndColumnMetaVariables.insert(variable);
}
+ std::sort(ddRowVariableIndices.begin(), ddRowVariableIndices.end());
for (auto const& variable : columnMetaVariables) {
DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
for (auto const& ddVariable : metaVariable.getDdVariables()) {
@@ -533,12 +541,14 @@ namespace storm {
}
rowAndColumnMetaVariables.insert(variable);
}
+ std::sort(ddColumnVariableIndices.begin(), ddColumnVariableIndices.end());
for (auto const& variable : groupMetaVariables) {
DdMetaVariable<DdType::CUDD> const& metaVariable = this->getDdManager()->getMetaVariable(variable);
for (auto const& ddVariable : metaVariable.getDdVariables()) {
ddGroupVariableIndices.push_back(ddVariable.getIndex());
}
}
+ std::sort(ddGroupVariableIndices.begin(), ddGroupVariableIndices.end());
// TODO: assert that the group variables are at the very top of the variable ordering?
@@ -705,7 +715,67 @@ namespace storm {
addToVectorRec(Cudd_T(dd), currentLevel + 1, maxLevel, currentOffset + odd.getElseOffset(), odd.getThenSuccessor(), ddVariableIndices, targetVector);
}
}
-
+
+ template<typename ValueType>
+ ADD Add<DdType::CUDD>::fromVector(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::set<storm::expressions::Variable> const& metaVariables) {
+ std::vector<uint_fast64_t> ddVariableIndices = getSortedVariableIndices(*ddManager, metaVariables);
+ uint_fast64_t offset = 0;
+ return ADD(ddManager->getCuddManager(), fromVectorRec(ddManager->getCuddManager().getManager(), offset, 0, ddVariableIndices.size(), values, odd, ddVariableIndices));
+ }
+
+ template<typename ValueType>
+ DdNode* Add<DdType::CUDD>::fromVectorRec(::DdManager* manager, uint_fast64_t& currentOffset, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices) {
+ if (currentLevel == maxLevel) {
+ // If we are in a terminal node of the ODD, we need to check whether the then-offset of the ODD is one
+ // (meaning the encoding is a valid one) or zero (meaning the encoding is not valid). Consequently, we
+ // need to copy the next value of the vector iff the then-offset is greater than zero.
+ if (odd.getThenOffset() > 0) {
+ return Cudd_addConst(manager, values[currentOffset++]);
+ } else {
+ return Cudd_ReadZero(manager);
+ }
+ } else {
+ // If the total offset is zero, we can just return the constant zero DD.
+ if (odd.getThenOffset() + odd.getElseOffset() == 0) {
+ return Cudd_ReadZero(manager);
+ }
+
+ // Determine the new else-successor.
+ DdNode* elseSuccessor = nullptr;
+ if (odd.getElseOffset() > 0) {
+ elseSuccessor = fromVectorRec(manager, currentOffset, currentLevel + 1, maxLevel, values, odd.getElseSuccessor(), ddVariableIndices);
+ } else {
+ elseSuccessor = Cudd_ReadZero(manager);
+ }
+ Cudd_Ref(elseSuccessor);
+
+ // Determine the new then-successor.
+ DdNode* thenSuccessor = nullptr;
+ if (odd.getThenOffset() > 0) {
+ thenSuccessor = fromVectorRec(manager, currentOffset, currentLevel + 1, maxLevel, values, odd.getThenSuccessor(), ddVariableIndices);
+ } else {
+ thenSuccessor = Cudd_ReadZero(manager);
+ }
+ Cudd_Ref(thenSuccessor);
+
+ // Create a node representing ITE(currentVar, thenSuccessor, elseSuccessor);
+ DdNode* result = Cudd_addIthVar(manager, static_cast<int>(ddVariableIndices[currentLevel]));
+ Cudd_Ref(result);
+ DdNode* newResult = Cudd_addIte(manager, result, thenSuccessor, elseSuccessor);
+ Cudd_Ref(newResult);
+
+ // Dispose of the intermediate results
+ Cudd_RecursiveDeref(manager, result);
+ Cudd_RecursiveDeref(manager, thenSuccessor);
+ Cudd_RecursiveDeref(manager, elseSuccessor);
+
+ // Before returning, we remove the protection imposed by the previous call to Cudd_Ref.
+ Cudd_Deref(newResult);
+
+ return newResult;
+ }
+ }
+
void Add<DdType::CUDD>::exportToDot(std::string const& filename) const {
if (filename.empty()) {
std::vector<ADD> cuddAddVector = { this->getCuddAdd() };
@@ -751,7 +821,7 @@ namespace storm {
DdForwardIterator<DdType::CUDD> Add<DdType::CUDD>::end(bool enumerateDontCareMetaVariables) const {
return DdForwardIterator<DdType::CUDD>(this->getDdManager(), nullptr, nullptr, 0, true, nullptr, enumerateDontCareMetaVariables);
}
-
+
std::ostream& operator<<(std::ostream& out, const Add<DdType::CUDD>& add) {
add.exportToDot();
return out;
diff --git a/src/storage/dd/CuddAdd.h b/src/storage/dd/CuddAdd.h
index 07cfcc20e..65de2dd10 100644
--- a/src/storage/dd/CuddAdd.h
+++ b/src/storage/dd/CuddAdd.h
@@ -27,6 +27,16 @@ namespace storm {
friend class Bdd<DdType::CUDD>;
friend class Odd<DdType::CUDD>;
+ /*!
+ * Creates an ADD from the given explicit vector.
+ *
+ * @param ddManager The manager responsible for this DD.
+ * @param values The vector that is to be represented by the ADD.
+ * @param odd An ODD that is used to do the translation.
+ * @param metaVariables The meta variables to use to encode the vector.
+ */
+ Add(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, std::vector<double> const& values, Odd<DdType::CUDD> const& odd, std::set<storm::expressions::Variable> const& metaVariables);
+
// Instantiate all copy/move constructors/assignments with the default implementation.
Add() = default;
Add(Add<DdType::CUDD> const& other) = default;
@@ -675,6 +685,33 @@ namespace storm {
template<typename ValueType>
void addToVectorRec(DdNode const* dd, uint_fast64_t currentLevel, uint_fast64_t maxLevel, uint_fast64_t currentOffset, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<ValueType>& targetVector) const;
+ /*!
+ * Builds an ADD representing the given vector.
+ *
+ * @param ddManager The manager responsible for the ADD.
+ * @param values The vector that is to be represented by the ADD.
+ * @param odd The ODD used for the translation.
+ * @param metaVariables The meta variables used for the translation.
+ * @return The resulting (CUDD) ADD.
+ */
+ template<typename ValueType>
+ static ADD fromVector(std::shared_ptr<DdManager<DdType::CUDD> const> ddManager, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::set<storm::expressions::Variable> const& metaVariables);
+
+ /*!
+ * Builds an ADD representing the given vector.
+ *
+ * @param manager The manager responsible for the ADD.
+ * @param currentOffset The current offset in the vector.
+ * @param currentLevel The current level in the DD.
+ * @param maxLevel The maximal level in the DD.
+ * @param values The vector that is to be represented by the ADD.
+ * @param odd The ODD used for the translation.
+ * @param ddVariableIndices The (sorted) list of DD variable indices to use.
+ * @return The resulting (CUDD) ADD node.
+ */
+ template<typename ValueType>
+ static DdNode* fromVectorRec(::DdManager* manager, uint_fast64_t& currentOffset, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<ValueType> const& values, Odd<DdType::CUDD> const& odd, std::vector<uint_fast64_t> const& ddVariableIndices);
+
// The ADD created by CUDD.
ADD cuddAdd;
};
diff --git a/src/storage/dd/CuddDd.cpp b/src/storage/dd/CuddDd.cpp
index 43e2afd3f..d56589824 100644
--- a/src/storage/dd/CuddDd.cpp
+++ b/src/storage/dd/CuddDd.cpp
@@ -50,9 +50,13 @@ namespace storm {
}
std::vector<uint_fast64_t> Dd<DdType::CUDD>::getSortedVariableIndices() const {
+ return getSortedVariableIndices(*this->getDdManager(), this->getContainedMetaVariables());
+ }
+
+ std::vector<uint_fast64_t> Dd<DdType::CUDD>::getSortedVariableIndices(DdManager<DdType::CUDD> const& manager, std::set<storm::expressions::Variable> const& metaVariables) {
std::vector<uint_fast64_t> ddVariableIndices;
- for (auto const& metaVariableName : this->getContainedMetaVariables()) {
- auto const& metaVariable = this->getDdManager()->getMetaVariable(metaVariableName);
+ for (auto const& metaVariableName : metaVariables) {
+ auto const& metaVariable = manager.getMetaVariable(metaVariableName);
for (auto const& ddVariable : metaVariable.getDdVariables()) {
ddVariableIndices.push_back(ddVariable.getIndex());
}
diff --git a/src/storage/dd/CuddDd.h b/src/storage/dd/CuddDd.h
index fdeb085ef..b80bfc671 100644
--- a/src/storage/dd/CuddDd.h
+++ b/src/storage/dd/CuddDd.h
@@ -105,6 +105,8 @@ namespace storm {
*/
std::shared_ptr<DdManager<DdType::CUDD> const> getDdManager() const;
+ protected:
+
/*!
* Retrieves the (sorted) list of the variable indices of DD variables contained in this DD.
*
@@ -112,7 +114,14 @@ namespace storm {
*/
std::vector<uint_fast64_t> getSortedVariableIndices() const;
- protected:
+ /*!
+ * Retrieves the (sorted) list of the variable indices of the DD variables given by the meta variable set.
+ *
+ * @param manager The manager responsible for the DD.
+ * @param metaVariable The set of meta variables for which to retrieve the index list.
+ * @return The sorted list of variable indices.
+ */
+ static std::vector<uint_fast64_t> getSortedVariableIndices(DdManager<DdType::CUDD> const& manager, std::set<storm::expressions::Variable> const& metaVariables);
/*!
* Retrieves the set of all meta variables contained in the DD.
diff --git a/src/storage/dd/CuddDdManager.cpp b/src/storage/dd/CuddDdManager.cpp
index 63a84a58b..d7adf1a6c 100644
--- a/src/storage/dd/CuddDdManager.cpp
+++ b/src/storage/dd/CuddDdManager.cpp
@@ -278,5 +278,9 @@ namespace storm {
void DdManager<DdType::CUDD>::triggerReordering() {
this->getCuddManager().ReduceHeap(this->reorderingTechnique, 0);
}
+
+ std::shared_ptr<DdManager<DdType::CUDD> const> DdManager<DdType::CUDD>::asSharedPointer() const {
+ return this->shared_from_this();
+ }
}
}
\ No newline at end of file
diff --git a/src/storage/dd/CuddDdManager.h b/src/storage/dd/CuddDdManager.h
index 07c11eb19..522be1fba 100644
--- a/src/storage/dd/CuddDdManager.h
+++ b/src/storage/dd/CuddDdManager.h
@@ -165,6 +165,13 @@ namespace storm {
*/
DdMetaVariable<DdType::CUDD> const& getMetaVariable(storm::expressions::Variable const& variable) const;
+ /*!
+ * Retrieves the manager as a shared pointer.
+ *
+ * @return A shared pointer to the manager.
+ */
+ std::shared_ptr<DdManager<DdType::CUDD> const> asSharedPointer() const;
+
private:
/*!
* Retrieves a list of names of the DD variables in the order of their index.
diff --git a/src/storage/dd/CuddOdd.cpp b/src/storage/dd/CuddOdd.cpp
index 40805d9e9..31bc42dd6 100644
--- a/src/storage/dd/CuddOdd.cpp
+++ b/src/storage/dd/CuddOdd.cpp
@@ -91,6 +91,10 @@ namespace storm {
}
}
+ bool Odd<DdType::CUDD>::isTerminalNode() const {
+ return this->elseNode == nullptr && this->thenNode == nullptr;
+ }
+
std::shared_ptr<Odd<DdType::CUDD>> Odd<DdType::CUDD>::buildOddFromAddRec(DdNode* dd, Cudd const& manager, uint_fast64_t currentLevel, uint_fast64_t maxLevel, std::vector<uint_fast64_t> const& ddVariableIndices, std::vector<std::unordered_map<DdNode*, std::shared_ptr<Odd<DdType::CUDD>>>>& uniqueTableForLevels) {
// Check whether the ODD for this node has already been computed (for this level) and if so, return this instead.
auto const& iterator = uniqueTableForLevels[currentLevel].find(dd);
diff --git a/src/storage/dd/CuddOdd.h b/src/storage/dd/CuddOdd.h
index 2d64275a3..9269e4ab1 100644
--- a/src/storage/dd/CuddOdd.h
+++ b/src/storage/dd/CuddOdd.h
@@ -96,6 +96,13 @@ namespace storm {
*/
uint_fast64_t getNodeCount() const;
+ /*!
+ * Checks whether the given ODD node is a terminal node, i.e. has no successors.
+ *
+ * @return True iff the node is terminal.
+ */
+ bool isTerminalNode() const;
+
private:
// Declare a hash functor that is used for the unique tables in the construction process.
class HashFunctor {
diff --git a/src/utility/cli.h b/src/utility/cli.h
index 202786ba2..780201cf3 100644
--- a/src/utility/cli.h
+++ b/src/utility/cli.h
@@ -501,7 +501,7 @@ namespace storm {
std::shared_ptr<storm::models::symbolic::Dtmc<DdType>> dtmc = model->template as<storm::models::symbolic::Dtmc<DdType>>();
storm::modelchecker::HybridDtmcPrctlModelChecker<DdType, double> modelchecker(*dtmc);
if (modelchecker.canHandle(*formula.get())) {
- modelchecker.check(*formula.get());
+ result = modelchecker.check(*formula.get());
}
} else {
STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This functionality is not yet implemented.");
diff --git a/src/utility/vector.h b/src/utility/vector.h
index bfc60b3a6..b5567ae5b 100644
--- a/src/utility/vector.h
+++ b/src/utility/vector.h
@@ -130,67 +130,78 @@ namespace storm {
}
/*!
- * Applies the given operation pointwise on the two given vectors and writes the result into the first
- * vector.
+ * Applies the given operation pointwise on the two given vectors and writes the result to the third vector.
+ * To botain an in-place operation, the third vector may be equal to any of the other two vectors.
*
- * @param target The first operand and target vector.
+ * @param firstOperand The first operand.
* @param secondOperand The second operand.
+ * @param target The target vector.
*/
template<class T>
- void applyPointwiseInPlace(std::vector<T>& target, std::vector<T> const& secondOperand, std::function<T (T const&, T const&)> function) {
-#ifdef DEBUG
- if (target.size() != summand.size()) {
- throw storm::exceptions::InvalidArgumentException() << "Invalid call to storm::utility::vector::applyPointwiseInPlace: operand lengths mismatch.";
- }
-#endif
+ void applyPointwise(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand, std::vector<T>& target, std::function<T (T const&, T const&)> function) {
#ifdef STORM_HAVE_INTELTBB
tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
[&](tbb::blocked_range<uint_fast64_t> const& range) {
- std::transform(target.begin() + range.begin(), target.begin() + range.end(), secondOperand.begin() + range.begin(), target.begin() + range.begin(), function);
+ std::transform(firstOperand.begin() + range.begin(), firstOperand.begin() + range.end(), secondOperand.begin() + range.begin(), target.begin() + range.begin(), function);
});
#else
- std::transform(target.begin(), target.end(), secondOperand.begin(), target.begin(), function);
+ std::transform(firstOperand.begin(), firstOperand.end(), secondOperand.begin(), target.begin(), function);
#endif
}
/*!
* Applies the given function pointwise on the given vector.
*
- * @param target The vector to which to apply the function.
+ * @param operand The vector to which to apply the function.
+ * @param target The target vector.
* @param function The function to apply.
*/
template<class T>
- void applyPointwiseInPlace(std::vector<T>& target, std::function<T (T const&)> const& function) {
+ void applyPointwise(std::vector<T> const& operand, std::vector<T>& target, std::function<T (T const&)> const& function) {
#ifdef STORM_HAVE_INTELTBB
tbb::parallel_for(tbb::blocked_range<uint_fast64_t>(0, target.size()),
[&](tbb::blocked_range<uint_fast64_t> const& range) {
- std::transform(target.begin() + range.begin(), target.begin() + range.end(), target.begin() + range.begin(), function);
+ std::transform(operand.begin() + range.begin(), operand.begin() + range.end(), target.begin() + range.begin(), function);
});
#else
- std::transform(target.begin(), target.end(), target.begin(), function);
+ std::transform(operand.begin(), operand.end(), target.begin(), function);
#endif
}
/*!
- * Adds the two given vectors and writes the result into the first operand.
+ * Adds the two given vectors and writes the result to the target vector.
*
- * @param target The first summand and target vector.
- * @param summand The second summand.
+ * @param firstOperand The first operand.
+ * @param secondOperand The second operand
+ * @param target The target vector.
*/
template<class T>
- void addVectorsInPlace(std::vector<T>& target, std::vector<T> const& summand) {
- applyPointwiseInPlace<T>(target, summand, std::plus<T>());
+ void addVectors(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand, std::vector<T>& target) {
+ applyPointwise<T>(firstOperand, secondOperand, target, std::plus<T>());
}
/*!
- * Subtracts the two given vectors and writes the result into the first operand.
+ * Subtracts the two given vectors and writes the result to the target vector.
*
- * @param target The first summand and target vector.
- * @param summand The second summand.
+ * @param firstOperand The first operand.
+ * @param secondOperand The second operand
+ * @param target The target vector.
+ */
+ template<class T>
+ void subtractVectors(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand, std::vector<T>& target) {
+ applyPointwise<T>(firstOperand, secondOperand, target, std::minus<T>());
+ }
+
+ /*!
+ * Multiplies the two given vectors (pointwise) and writes the result to the target vector.
+ *
+ * @param firstOperand The first operand.
+ * @param secondOperand The second operand
+ * @param target The target vector.
*/
template<class T>
- void subtractVectorsInPlace(std::vector<T>& target, std::vector<T> const& summand) {
- applyPointwiseInPlace<T>(target, summand, std::minus<T>());
+ void multiplyVectorsPointwise(std::vector<T> const& firstOperand, std::vector<T> const& secondOperand, std::vector<T>& target) {
+ applyPointwise<T>(firstOperand, secondOperand, target, std::multiplies<T>());
}
/*!
@@ -201,7 +212,7 @@ namespace storm {
*/
template<class T>
void scaleVectorInPlace(std::vector<T>& target, T const& factor) {
- applyPointwiseInPlace<T>(target, [&] (T const& argument) { return argument * factor; });
+ applyPointwise<T>(target, target, [&] (T const& argument) { return argument * factor; });
}
/*!
diff --git a/test/functional/storage/CuddDdTest.cpp b/test/functional/storage/CuddDdTest.cpp
index 71efcadaf..6255e912c 100644
--- a/test/functional/storage/CuddDdTest.cpp
+++ b/test/functional/storage/CuddDdTest.cpp
@@ -376,6 +376,9 @@ TEST(CuddDd, BddOddTest) {
EXPECT_TRUE(i+1 == ddAsVector[i]);
}
+ storm::dd::Add<storm::dd::DdType::CUDD> vectorAdd(manager, ddAsVector, odd, {x.first});
+ EXPECT_EQ(dd, vectorAdd);
+
// Create a non-trivial matrix.
dd = manager->getIdentity(x.first).equals(manager->getIdentity(x.second)) * manager->getRange(x.first).toAdd();
dd += manager->getEncoding(x.first, 1).toAdd() * manager->getRange(x.second).toAdd() + manager->getEncoding(x.second, 1).toAdd() * manager->getRange(x.first).toAdd();
diff --git a/test/functional/storage/SparseMatrixTest.cpp b/test/functional/storage/SparseMatrixTest.cpp
index b677f6f3c..20716916a 100644
--- a/test/functional/storage/SparseMatrixTest.cpp
+++ b/test/functional/storage/SparseMatrixTest.cpp
@@ -447,7 +447,7 @@ TEST(SparseMatrix, JacobiDecomposition) {
ASSERT_NO_THROW(matrix = matrixBuilder.build());
ASSERT_NO_THROW(matrix.getJacobiDecomposition());
- std::pair<storm::storage::SparseMatrix<double>, storm::storage::SparseMatrix<double>> jacobiDecomposition = matrix.getJacobiDecomposition();
+ std::pair<storm::storage::SparseMatrix<double>, std::vector<double>> jacobiDecomposition = matrix.getJacobiDecomposition();
storm::storage::SparseMatrixBuilder<double> luBuilder(4, 4, 3);
ASSERT_NO_THROW(luBuilder.addNextValue(0, 1, 1.2));
@@ -456,14 +456,7 @@ TEST(SparseMatrix, JacobiDecomposition) {
storm::storage::SparseMatrix<double> lu;
ASSERT_NO_THROW(lu = luBuilder.build());
- storm::storage::SparseMatrixBuilder<double> dinvBuilder(4, 4, 4);
- ASSERT_NO_THROW(dinvBuilder.addNextValue(0, 0, 1 / 1.1));
- ASSERT_NO_THROW(dinvBuilder.addNextValue(1, 1, 1 / 0.5));
- ASSERT_NO_THROW(dinvBuilder.addNextValue(2, 2, 1 / 0.99));
- ASSERT_NO_THROW(dinvBuilder.addNextValue(3, 3, 1 / 0.11));
- storm::storage::SparseMatrix<double> dinv;
- ASSERT_NO_THROW(dinv = dinvBuilder.build());
-
+ std::vector<double> dinv = {1/1.1, 1/0.5, 1/0.99, 1/0.11};
ASSERT_TRUE(lu == jacobiDecomposition.first);
ASSERT_TRUE(dinv == jacobiDecomposition.second);
}