Moved model checking functionality for MDPs for general superclass such that specialized model checkers only need to implement certain operations. Fixed tests.

src/modelchecker/DtmcPrctlModelChecker.h

@@ -115,7 +115,7 @@ public:
 		storm::utility::setVectorValues(result, *rightStates, storm::utility::constGetOne<Type>());
 
 		// Perform the matrix vector multiplication as often as required by the formula bound.
-		this->performMatrixVectorMultiplication(tmpMatrix, &result, nullptr, formula.getBound());
+		this->performMatrixVectorMultiplication(tmpMatrix, *result, nullptr, formula.getBound());
 
 		// Return result.
 		return result;
@@ -139,7 +139,7 @@ public:
 		delete nextStates;
 
 		// Perform one single matrix-vector multiplication.
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), &result);
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result);
 
 		// Return result.
 		return result;
@@ -228,20 +228,20 @@ public:
 			// Initialize the x vector with 0.5 for each element. This is the initial guess for
 			// the iterative solvers. It should be safe as for all 'maybe' states we know that the
 			// probability is strictly larger than 0.
-			std::vector<Type>* x = new std::vector<Type>(maybeStatesSetBitCount, Type(0.5));
+			std::vector<Type> x(maybeStatesSetBitCount, Type(0.5));
 
 			// Prepare the right-hand side of the equation system. For entry i this corresponds to
 			// the accumulated probability of going from state i to some 'yes' state.
 			std::vector<Type> b(maybeStatesSetBitCount);
 			this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, statesWithProbability1, &b);
 
-			this->solveEquationSystem(*submatrix, &x, b);
+			this->solveEquationSystem(*submatrix, x, b);
 
 			// Delete the created submatrix.
 			delete submatrix;
 
 			// Set values of resulting vector according to result.
-			storm::utility::setVectorValues<Type>(result, maybeStates, *x);
+			storm::utility::setVectorValues<Type>(result, maybeStates, x);
 		} else if (qualitative) {
 			// If we only need a qualitative result, we can safely assume that the results will only be compared to
 			// bounds which are either 0 or 1. Setting the value to 0.5 is thus safe.
@@ -273,7 +273,7 @@ public:
 		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
 
 		// Perform the actual matrix-vector multiplication as long as the bound of the formula is met.
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), &result, nullptr, formula.getBound());
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, nullptr, formula.getBound());
 
 		// Return result.
 		return result;
@@ -306,7 +306,7 @@ public:
 
 		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), &result, totalRewardVector, formula.getBound());
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, totalRewardVector, formula.getBound());
 
 		// Delete temporary variables and return result.
 		delete totalRewardVector;
@@ -351,16 +351,16 @@ public:
 
 			// Initialize the x vector with 1 for each element. This is the initial guess for
 			// the iterative solvers.
-			std::vector<Type>* x = new std::vector<Type>(maybeStatesSetBitCount, storm::utility::constGetOne<Type>());
+			std::vector<Type> x(maybeStatesSetBitCount, storm::utility::constGetOne<Type>());
 
 			// Prepare the right-hand side of the equation system.
-			std::vector<Type>* b = new std::vector<Type>(maybeStatesSetBitCount);
+			std::vector<Type> b(maybeStatesSetBitCount);
 			if (this->getModel().hasTransitionRewards()) {
 				// If a transition-based reward model is available, we initialize the right-hand
 				// side to the vector resulting from summing the rows of the pointwise product
 				// of the transition probability matrix and the transition reward matrix.
 				std::vector<Type>* pointwiseProductRowSumVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
-				storm::utility::selectVectorValues(b, maybeStates, *pointwiseProductRowSumVector);
+				storm::utility::selectVectorValues(&b, maybeStates, *pointwiseProductRowSumVector);
 				delete pointwiseProductRowSumVector;
 
 				if (this->getModel().hasStateRewards()) {
@@ -368,27 +368,25 @@ public:
 					// as well. As the state reward vector contains entries not just for the states
 					// that we still consider (i.e. maybeStates), we need to extract these values
 					// first.
-					std::vector<Type>* subStateRewards = new std::vector<Type>(maybeStatesSetBitCount);
-					storm::utility::selectVectorValues(subStateRewards, maybeStates, *this->getModel().getStateRewardVector());
-					gmm::add(*subStateRewards, *b);
-					delete subStateRewards;
+					std::vector<Type> subStateRewards(maybeStatesSetBitCount);
+					storm::utility::selectVectorValues(&subStateRewards, maybeStates, *this->getModel().getStateRewardVector());
+					gmm::add(subStateRewards, b);
 				}
 			} else {
 				// If only a state-based reward model is  available, we take this vector as the
 				// right-hand side. As the state reward vector contains entries not just for the
 				// states that we still consider (i.e. maybeStates), we need to extract these values
 				// first.
-				storm::utility::selectVectorValues(b, maybeStates, *this->getModel().getStateRewardVector());
+				storm::utility::selectVectorValues(&b, maybeStates, *this->getModel().getStateRewardVector());
 			}
 
-			this->solveEquationSystem(*submatrix, &x, *b);
+			this->solveEquationSystem(*submatrix, x, b);
 
 			// Set values of resulting vector according to result.
-			storm::utility::setVectorValues<Type>(result, maybeStates, *x);
+			storm::utility::setVectorValues<Type>(result, maybeStates, x);
 
 			// Delete temporary matrix and right-hand side.
 			delete submatrix;
-			delete b;
 		}
 
 		// Set values of resulting vector that are known exactly.
@@ -401,9 +399,9 @@ public:
 	}
 
 private:
-	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>** vector, std::vector<Type>* summand = nullptr, uint_fast64_t repetitions = 1) const = 0;
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand = nullptr, uint_fast64_t repetitions = 1) const = 0;
 
-	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>** vector, std::vector<Type>& b) const = 0;
+	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type> const& b) const = 0;
 };
 
 } //namespace modelChecker

src/modelchecker/GmmxxDtmcPrctlModelChecker.h

@@ -96,25 +96,33 @@ public:
 
 private:
 
-	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>** vector, std::vector<Type>* summand, uint_fast64_t repetitions = 1) const {
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand, uint_fast64_t repetitions = 1) const {
 		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
 		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
 
 		// Now perform matrix-vector multiplication as long as we meet the bound.
 		std::vector<Type>* swap = nullptr;
-		std::vector<Type>* tmpResult = new std::vector<Type>(this->getModel().getNumberOfStates());
+		std::vector<Type>* currentVector = &vector;
+		std::vector<Type>* tmpVector = new std::vector<Type>(this->getModel().getNumberOfStates());
 		for (uint_fast64_t i = 0; i < repetitions; ++i) {
-			gmm::mult(*gmmxxMatrix, **vector, *tmpResult);
-			swap = tmpResult;
-			tmpResult = *vector;
-			*vector = swap;
+			gmm::mult(*gmmxxMatrix, *currentVector, *tmpVector);
+			swap = tmpVector;
+			tmpVector = currentVector;
+			currentVector = swap;
 
 			// If requested, add an offset to the current result vector.
 			if (summand != nullptr) {
-				gmm::add(*summand, **vector);
+				gmm::add(*summand, *currentVector);
 			}
 		}
-		delete tmpResult;
+
+		if (repetitions % 2 == 1) {
+			std::swap(vector, *currentVector);
+			delete currentVector;
+		} else {
+			delete tmpVector;
+		}
+
 		delete gmmxxMatrix;
 	}
 
@@ -128,7 +136,7 @@ private:
 	 * @return The solution of the system of linear equations in form of the elements of the vector
 	 * x.
 	 */
-	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>** vector, std::vector<Type>& b) const {
+	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type> const& b) const {
 		// Get the settings object to customize linear solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
@@ -155,13 +163,13 @@ private:
 		if (s->getString("lemethod") == "bicgstab") {
 			LOG4CPLUS_INFO(logger, "Using BiCGStab method.");
 			if (precond == "ilu") {
-				gmm::bicgstab(*A, **vector, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::bicgstab(*A, vector, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*A), iter);
 			} else if (precond == "diagonal") {
-				gmm::bicgstab(*A, **vector, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::bicgstab(*A, vector, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*A), iter);
 			} else if (precond == "ildlt") {
-				gmm::bicgstab(*A, **vector, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::bicgstab(*A, vector, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*A), iter);
 			} else if (precond == "none") {
-				gmm::bicgstab(*A, **vector, b, gmm::identity_matrix(), iter);
+				gmm::bicgstab(*A, vector, b, gmm::identity_matrix(), iter);
 			}
 
 			// Check if the solver converged and issue a warning otherwise.
@@ -173,13 +181,13 @@ private:
 		} else if (s->getString("lemethod") == "qmr") {
 			LOG4CPLUS_INFO(logger, "Using QMR method.");
 			if (precond == "ilu") {
-				gmm::qmr(*A, **vector, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::qmr(*A, vector, b, gmm::ilu_precond<gmm::csr_matrix<Type>>(*A), iter);
 			} else if (precond == "diagonal") {
-				gmm::qmr(*A, **vector, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::qmr(*A, vector, b, gmm::diagonal_precond<gmm::csr_matrix<Type>>(*A), iter);
 			} else if (precond == "ildlt") {
-				gmm::qmr(*A, **vector, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*A), iter);
+				gmm::qmr(*A, vector, b, gmm::ildlt_precond<gmm::csr_matrix<Type>>(*A), iter);
 			} else if (precond == "none") {
-				gmm::qmr(*A, **vector, b, gmm::identity_matrix(), iter);
+				gmm::qmr(*A, vector, b, gmm::identity_matrix(), iter);
 			}
 
 			// Check if the solver converged and issue a warning otherwise.
@@ -190,7 +198,7 @@ private:
 			}
 		} else if (s->getString("lemethod") == "jacobi") {
 			LOG4CPLUS_INFO(logger, "Using Jacobi method.");
-			solveLinearEquationSystemWithJacobi(*A, **vector, b);
+			solveLinearEquationSystemWithJacobi(*A, vector, b);
 		}
 
 		delete A;

src/modelchecker/GmmxxMdpPrctlModelChecker.h

@@ -43,339 +43,33 @@ public:
 
 	virtual ~GmmxxMdpPrctlModelChecker() { }
 
-	virtual std::vector<Type>* checkBoundedUntil(const storm::formula::BoundedUntil<Type>& formula, bool qualitative) const {
-		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
-		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
-		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
-
-		// Copy the matrix before we make any changes.
-		storm::storage::SparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionMatrix());
-
+private:
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand = nullptr, uint_fast64_t repetitions = 1) const {
 		// Get the starting row indices for the non-deterministic choices to reduce the resulting
 		// vector properly.
 		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
 
-		// Make all rows absorbing that violate both sub-formulas or satisfy the second sub-formula.
-		tmpMatrix.makeRowsAbsorbing(~(*leftStates | *rightStates) | *rightStates, *nondeterministicChoiceIndices);
-
 		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(tmpMatrix);
-
-		// Create the vector with which to multiply.
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
-		storm::utility::setVectorValues(result, *rightStates, storm::utility::constGetOne<Type>());
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
 
 		// Create vector for result of multiplication, which is reduced to the result vector after
 		// each multiplication.
-		std::vector<Type>* multiplyResult = new std::vector<Type>(this->getModel().getTransitionMatrix()->getRowCount(), 0);
+		std::vector<Type> multiplyResult(matrix.getRowCount());
 
 		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
-		for (uint_fast64_t i = 0; i < formula.getBound(); ++i) {
-			gmm::mult(*gmmxxMatrix, *result, *multiplyResult);
+		for (uint_fast64_t i = 0; i < repetitions; ++i) {
+			gmm::mult(*gmmxxMatrix, vector, multiplyResult);
 			if (this->minimumOperatorStack.top()) {
-				storm::utility::reduceVectorMin(*multiplyResult, result, *nondeterministicChoiceIndices);
+				storm::utility::reduceVectorMin(multiplyResult, &vector, *nondeterministicChoiceIndices);
 			} else {
-				storm::utility::reduceVectorMax(*multiplyResult, result, *nondeterministicChoiceIndices);
+				storm::utility::reduceVectorMax(multiplyResult, &vector, *nondeterministicChoiceIndices);
 			}
 		}
-		delete multiplyResult;
 
 		// Delete intermediate results and return result.
 		delete gmmxxMatrix;
-		delete leftStates;
-		delete rightStates;
-		return result;
-	}
-
-	virtual std::vector<Type>* checkNext(const storm::formula::Next<Type>& formula, bool qualitative) const {
-		// First, we need to compute the states that satisfy the sub-formula of the next-formula.
-		storm::storage::BitVector* nextStates = formula.getChild().check(*this);
-
-		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*this->getModel().getTransitionMatrix());
-
-		// Create the vector with which to multiply and initialize it correctly.
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
-		storm::utility::setVectorValues(result, *nextStates, storm::utility::constGetOne<Type>());
-
-		// Delete obsolete sub-result.
-		delete nextStates;
-
-		// Create resulting vector.
-		std::vector<Type>* temporaryResult = new std::vector<Type>(this->getModel().getTransitionMatrix()->getRowCount());
-
-		// Perform the actual computation, namely matrix-vector multiplication.
-		gmm::mult(*gmmxxMatrix, *result, *temporaryResult);
-
-		// Get the starting row indices for the non-deterministic choices to reduce the resulting
-		// vector properly.
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
-
-		if (this->minimumOperatorStack.top()) {
-			storm::utility::reduceVectorMin(*temporaryResult, result, *nondeterministicChoiceIndices);
-		} else {
-			storm::utility::reduceVectorMax(*temporaryResult, result, *nondeterministicChoiceIndices);
-		}
-
-		// Delete temporary matrix plus temporary result and return result.
-		delete gmmxxMatrix;
-		delete temporaryResult;
-		return result;
-	}
-
-	virtual std::vector<Type>* checkUntil(const storm::formula::Until<Type>& formula, bool qualitative) const {
-		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
-		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
-		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
-
-		// Then, we need to identify the states which have to be taken out of the matrix, i.e.
-		// all states that have probability 0 and 1 of satisfying the until-formula.
-		storm::storage::BitVector statesWithProbability0(this->getModel().getNumberOfStates());
-		storm::storage::BitVector statesWithProbability1(this->getModel().getNumberOfStates());
-		if (this->minimumOperatorStack.top()) {
-			storm::utility::GraphAnalyzer::performProb01Min(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
-		} else {
-			storm::utility::GraphAnalyzer::performProb01Max(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
-		}
-
-		// Delete sub-results that are obsolete now.
-		delete leftStates;
-		delete rightStates;
-
-		LOG4CPLUS_INFO(logger, "Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states.");
-		LOG4CPLUS_INFO(logger, "Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states.");
-		storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
-		LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
-
-		// Create resulting vector.
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
-
-		// Only try to solve system if there are states for which the probability is unknown.
-		uint_fast64_t maybeStatesSetBitCount = maybeStates.getNumberOfSetBits();
-		if (maybeStatesSetBitCount > 0) {
-			// First, we can eliminate the rows and columns from the original transition probability matrix for states
-			// whose probabilities are already known.
-			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates, *this->getModel().getNondeterministicChoiceIndices());
-
-			// Get the "new" nondeterministic choice indices for the submatrix.
-			std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
-
-			// Transform the submatrix to the gmm++ format to use its capabilities.
-			gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*submatrix);
-
-			// Create vector for results for maybe states.
-			std::vector<Type>* x = new std::vector<Type>(maybeStatesSetBitCount);
-
-			// Prepare the right-hand side of the equation system. For entry i this corresponds to
-			// the accumulated probability of going from state i to some 'yes' state.
-			std::vector<Type> b(submatrix->getRowCount());
-			this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, *this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, &b);
-			delete submatrix;
-
-			// Solve the corresponding system of equations.
-			this->solveEquationSystem(*gmmxxMatrix, x, b, *subNondeterministicChoiceIndices);
-
-			// Set values of resulting vector according to result.
-			storm::utility::setVectorValues<Type>(result, maybeStates, *x);
-
-			// Delete temporary matrix and vector.
-			delete gmmxxMatrix;
-			delete x;
-		}
-
-		// Set values of resulting vector that are known exactly.
-		storm::utility::setVectorValues<Type>(result, statesWithProbability0, storm::utility::constGetZero<Type>());
-		storm::utility::setVectorValues<Type>(result, statesWithProbability1, storm::utility::constGetOne<Type>());
-
-		return result;
-	}
-
-	virtual std::vector<Type>* checkInstantaneousReward(const storm::formula::InstantaneousReward<Type>& formula, bool qualitative) const {
-		// Only compute the result if the model has a state-based reward model.
-		if (!this->getModel().hasStateRewards()) {
-			LOG4CPLUS_ERROR(logger, "Missing (state-based) reward model for formula.");
-			throw storm::exceptions::InvalidPropertyException() << "Missing (state-based) reward model for formula.";
-		}
-
-		// Get the starting row indices for the non-deterministic choices to reduce the resulting
-		// vector properly.
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
-
-		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*this->getModel().getTransitionMatrix());
-
-		// Initialize result to state rewards of the model.
-		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
-
-		// Create vector for result of multiplication, which is reduced to the result vector after
-		// each multiplication.
-		std::vector<Type>* multiplyResult = new std::vector<Type>(this->getModel().getTransitionMatrix()->getRowCount(), 0);
-
-		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
-		for (uint_fast64_t i = 0; i < formula.getBound(); ++i) {
-			gmm::mult(*gmmxxMatrix, *result, *multiplyResult);
-			if (this->minimumOperatorStack.top()) {
-				storm::utility::reduceVectorMin(*multiplyResult, result, *nondeterministicChoiceIndices);
-			} else {
-				storm::utility::reduceVectorMax(*multiplyResult, result, *nondeterministicChoiceIndices);
-			}
-		}
-		delete multiplyResult;
-
-		// Delete temporary variables and return result.
-		delete gmmxxMatrix;
-		return result;
-	}
-
-	virtual std::vector<Type>* checkCumulativeReward(const storm::formula::CumulativeReward<Type>& formula, bool qualitative) const {
-		// Only compute the result if the model has at least one reward model.
-		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
-			LOG4CPLUS_ERROR(logger, "Missing reward model for formula.");
-			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
-		}
-
-		// Get the starting row indices for the non-deterministic choices to reduce the resulting
-		// vector properly.
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
-
-		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
-		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*this->getModel().getTransitionMatrix());
-
-		// Compute the reward vector to add in each step based on the available reward models.
-		std::vector<Type>* totalRewardVector = nullptr;
-		if (this->getModel().hasTransitionRewards()) {
-			totalRewardVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
-			if (this->getModel().hasStateRewards()) {
-				gmm::add(*this->getModel().getStateRewardVector(), *totalRewardVector);
-			}
-		} else {
-			totalRewardVector = new std::vector<Type>(*this->getModel().getStateRewardVector());
-		}
-
-		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
-
-		// Create vector for result of multiplication, which is reduced to the result vector after
-		// each multiplication.
-		std::vector<Type>* multiplyResult = new std::vector<Type>(this->getModel().getTransitionMatrix()->getRowCount(), 0);
-
-		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
-		for (uint_fast64_t i = 0; i < formula.getBound(); ++i) {
-			gmm::mult(*gmmxxMatrix, *result, *multiplyResult);
-			if (this->minimumOperatorStack.top()) {
-				storm::utility::reduceVectorMin(*multiplyResult, result, *nondeterministicChoiceIndices);
-			} else {
-				storm::utility::reduceVectorMax(*multiplyResult, result, *nondeterministicChoiceIndices);
-			}
-
-			// Add the reward vector to the result.
-			gmm::add(*totalRewardVector, *result);
-		}
-		delete multiplyResult;
-
-		// Delete temporary variables and return result.
-		delete gmmxxMatrix;
-		delete totalRewardVector;
-		return result;
 	}
 
-	virtual std::vector<Type>* checkReachabilityReward(const storm::formula::ReachabilityReward<Type>& formula, bool qualitative) const {
-		// Only compute the result if the model has at least one reward model.
-		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
-			LOG4CPLUS_ERROR(logger, "Missing reward model for formula. Skipping formula");
-			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
-		}
-
-		// Determine the states for which the target predicate holds.
-		storm::storage::BitVector* targetStates = formula.getChild().check(*this);
-
-		// Determine which states have a reward of infinity by definition.
-		storm::storage::BitVector infinityStates(this->getModel().getNumberOfStates());
-		storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true);
-		if (this->minimumOperatorStack.top()) {
-			storm::utility::GraphAnalyzer::performProb1A(this->getModel(), trueStates, *targetStates, &infinityStates);
-		} else {
-			storm::utility::GraphAnalyzer::performProb1E(this->getModel(), trueStates, *targetStates, &infinityStates);
-		}
-		infinityStates.complement();
-
-		LOG4CPLUS_INFO(logger, "Found " << infinityStates.getNumberOfSetBits() << " 'infinity' states.");
-		LOG4CPLUS_INFO(logger, "Found " << targetStates->getNumberOfSetBits() << " 'target' states.");
-		storm::storage::BitVector maybeStates = ~(*targetStates) & ~infinityStates;
-		LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
-
-		// Create resulting vector.
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
-
-		// Get the starting row indices for the non-deterministic choices to reduce the resulting
-		// vector properly.
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
-
-		// Check whether there are states for which we have to compute the result.
-		const int maybeStatesSetBitCount = maybeStates.getNumberOfSetBits();
-		if (maybeStatesSetBitCount > 0) {
-			// First, we can eliminate the rows and columns from the original transition probability matrix for states
-			// whose probabilities are already known.
-			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates, *this->getModel().getNondeterministicChoiceIndices());
-
-			// Get the "new" nondeterministic choice indices for the submatrix.
-			std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
-
-			// Transform the submatrix to the gmm++ format to use its capabilities.
-			gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*submatrix);
-
-			// Create vector for results for maybe states.
-			std::vector<Type>* x = new std::vector<Type>(maybeStatesSetBitCount);
-
-			// Prepare the right-hand side of the equation system. For entry i this corresponds to
-			// the accumulated probability of going from state i to some 'yes' state.
-			std::vector<Type> b(submatrix->getRowCount());
-			delete submatrix;
-
-			if (this->getModel().hasTransitionRewards()) {
-				// If a transition-based reward model is available, we initialize the right-hand
-				// side to the vector resulting from summing the rows of the pointwise product
-				// of the transition probability matrix and the transition reward matrix.
-				std::vector<Type>* pointwiseProductRowSumVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
-				storm::utility::selectVectorValues(&b, maybeStates, *nondeterministicChoiceIndices, *pointwiseProductRowSumVector);
-				delete pointwiseProductRowSumVector;
-
-				if (this->getModel().hasStateRewards()) {
-					// If a state-based reward model is also available, we need to add this vector
-					// as well. As the state reward vector contains entries not just for the states
-					// that we still consider (i.e. maybeStates), we need to extract these values
-					// first.
-					std::vector<Type>* subStateRewards = new std::vector<Type>(b.size());
-					storm::utility::selectVectorValuesRepeatedly(subStateRewards, maybeStates, *nondeterministicChoiceIndices, *this->getModel().getStateRewardVector());
-					gmm::add(*subStateRewards, b);
-					delete subStateRewards;
-				}
-			} else {
-				// If only a state-based reward model is  available, we take this vector as the
-				// right-hand side. As the state reward vector contains entries not just for the
-				// states that we still consider (i.e. maybeStates), we need to extract these values
-				// first.
-				storm::utility::selectVectorValuesRepeatedly(&b, maybeStates, *nondeterministicChoiceIndices, *this->getModel().getStateRewardVector());
-			}
-
-			// Solve the corresponding system of equations.
-			this->solveEquationSystem(*gmmxxMatrix, x, b, *subNondeterministicChoiceIndices);
-
-			// Set values of resulting vector according to result.
-			storm::utility::setVectorValues<Type>(result, maybeStates, *x);
-			delete x;
-			delete gmmxxMatrix;
-		}
-
-		// Set values of resulting vector that are known exactly.
-		storm::utility::setVectorValues(result, *targetStates, storm::utility::constGetZero<Type>());
-		storm::utility::setVectorValues(result, infinityStates, storm::utility::constGetInfinity<Type>());
-
-		// Delete temporary storages and return result.
-		delete targetStates;
-		return result;
-	}
-
-private:
 	/*!
 	 * Solves the given equation system under the given parameters using the power method.
 	 *
@@ -386,7 +80,7 @@ private:
 	 * @return The solution of the system of linear equations in form of the elements of the vector
 	 * x.
 	 */
-	void solveEquationSystem(gmm::csr_matrix<Type> const& A, std::vector<Type>* x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
+	void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& x, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const {
 		// Get the settings object to customize solving.
 		storm::settings::Settings* s = storm::settings::instance();
 
@@ -395,43 +89,48 @@ private:
 		unsigned maxIterations = s->get<unsigned>("maxiter");
 		bool relative = s->get<bool>("relative");
 
+		// Transform the transition probability matrix to the gmm++ format to use its arithmetic.
+		gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(matrix);
+
 		// Set up the environment for the power method.
-		std::vector<Type>* temporaryResult = new std::vector<Type>(b.size());
-		std::vector<Type>* newX = new std::vector<Type>(x->size());
+		std::vector<Type> multiplyResult(matrix.getRowCount());
+		std::vector<Type>* currentX = &x;
+		std::vector<Type>* newX = new std::vector<Type>(x.size());
 		std::vector<Type>* swap = nullptr;
-		bool converged = false;
 		uint_fast64_t iterations = 0;
+		bool converged = false;
 
 		// Proceed with the iterations as long as the method did not converge or reach the
 		// user-specified maximum number of iterations.
 		while (!converged && iterations < maxIterations) {
 			// Compute x' = A*x + b.
-			gmm::mult(A, *x, *temporaryResult);
-			gmm::add(b, *temporaryResult);
+			gmm::mult(*gmmxxMatrix, *currentX, multiplyResult);
+			gmm::add(b, multiplyResult);
 
 			// Reduce the vector x' by applying min/max for all non-deterministic choices.
 			if (this->minimumOperatorStack.top()) {
-				storm::utility::reduceVectorMin(*temporaryResult, newX, nondeterministicChoiceIndices);
+				storm::utility::reduceVectorMin(multiplyResult, newX, nondeterministicChoiceIndices);
 			} else {
-				storm::utility::reduceVectorMax(*temporaryResult, newX, nondeterministicChoiceIndices);
+				storm::utility::reduceVectorMax(multiplyResult, newX, nondeterministicChoiceIndices);
 			}
 
 			// Determine whether the method converged.
-			converged = storm::utility::equalModuloPrecision(*x, *newX, precision, relative);
+			converged = storm::utility::equalModuloPrecision(*currentX, *newX, precision, relative);
 
 			// Update environment variables.
-			swap = x;
-			x = newX;
+			swap = currentX;
+			currentX = newX;
 			newX = swap;
 			++iterations;
 		}
 
 		if (iterations % 2 == 1) {
-			delete x;
+			std::swap(x, *currentX);
+			delete currentX;
 		} else {
 			delete newX;
 		}
-		delete temporaryResult;
+		delete gmmxxMatrix;
 
 		// Check if the solver converged and issue a warning otherwise.
 		if (converged) {
@@ -440,22 +139,6 @@ private:
 			LOG4CPLUS_WARN(logger, "Iterative solver did not converge.");
 		}
 	}
-
-	std::shared_ptr<std::vector<uint_fast64_t>> computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector constraint) const {
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
-		std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices(new std::vector<uint_fast64_t>(constraint.getNumberOfSetBits() + 1));
-		uint_fast64_t currentRowCount = 0;
-		uint_fast64_t currentIndexCount = 1;
-		(*subNondeterministicChoiceIndices)[0] = 0;
-		for (auto index : constraint) {
-			(*subNondeterministicChoiceIndices)[currentIndexCount] = currentRowCount + (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index];
-			currentRowCount += (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index];
-			++currentIndexCount;
-		}
-		(*subNondeterministicChoiceIndices)[constraint.getNumberOfSetBits()] = currentRowCount;
-
-		return subNondeterministicChoiceIndices;
-	}
 };
 
 } //namespace modelChecker

src/modelchecker/MdpPrctlModelChecker.h

@@ -189,7 +189,28 @@ public:
 	 * @param formula The Bounded Until path formula to check
 	 * @returns for each state the probability that the path formula holds.
 	 */
-	virtual std::vector<Type>* checkBoundedUntil(const storm::formula::BoundedUntil<Type>& formula, bool qualitative) const = 0;
+	virtual std::vector<Type>* checkBoundedUntil(const storm::formula::BoundedUntil<Type>& formula, bool qualitative) const {
+		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
+		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
+		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
+
+		// Copy the matrix before we make any changes.
+		storm::storage::SparseMatrix<Type> tmpMatrix(*this->getModel().getTransitionMatrix());
+
+		// Make all rows absorbing that violate both sub-formulas or satisfy the second sub-formula.
+		tmpMatrix.makeRowsAbsorbing(~(*leftStates | *rightStates) | *rightStates, *this->getModel().getNondeterministicChoiceIndices());
+
+		// Create the vector with which to multiply.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		storm::utility::setVectorValues(result, *rightStates, storm::utility::constGetOne<Type>());
+
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, nullptr, formula.getBound());
+
+		// Delete intermediate results and return result.
+		delete leftStates;
+		delete rightStates;
+		return result;
+	}
 
 	/*!
 	 * The check method for a path formula with a Next operator node as root in its formula tree
@@ -197,7 +218,22 @@ public:
 	 * @param formula The Next path formula to check
 	 * @returns for each state the probability that the path formula holds.
 	 */
-	virtual std::vector<Type>* checkNext(const storm::formula::Next<Type>& formula, bool qualitative) const = 0;
+	virtual std::vector<Type>* checkNext(const storm::formula::Next<Type>& formula, bool qualitative) const {
+		// First, we need to compute the states that satisfy the sub-formula of the next-formula.
+		storm::storage::BitVector* nextStates = formula.getChild().check(*this);
+
+		// Create the vector with which to multiply and initialize it correctly.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		storm::utility::setVectorValues(result, *nextStates, storm::utility::constGetOne<Type>());
+
+		// Delete obsolete sub-result.
+		delete nextStates;
+
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result);
+
+		// Return result.
+		return result;
+	}
 
 	/*!
 	 * The check method for a path formula with a Bounded Eventually operator node as root in its
@@ -246,7 +282,67 @@ public:
 	 * @param formula The Until path formula to check
 	 * @returns for each state the probability that the path formula holds.
 	 */
-	virtual std::vector<Type>* checkUntil(const storm::formula::Until<Type>& formula, bool qualitative) const = 0;
+	virtual std::vector<Type>* checkUntil(const storm::formula::Until<Type>& formula, bool qualitative) const {
+		// First, we need to compute the states that satisfy the sub-formulas of the until-formula.
+		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
+		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
+
+		// Then, we need to identify the states which have to be taken out of the matrix, i.e.
+		// all states that have probability 0 and 1 of satisfying the until-formula.
+		storm::storage::BitVector statesWithProbability0(this->getModel().getNumberOfStates());
+		storm::storage::BitVector statesWithProbability1(this->getModel().getNumberOfStates());
+		if (this->minimumOperatorStack.top()) {
+			storm::utility::GraphAnalyzer::performProb01Min(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
+		} else {
+			storm::utility::GraphAnalyzer::performProb01Max(this->getModel(), *leftStates, *rightStates, &statesWithProbability0, &statesWithProbability1);
+		}
+
+		// Delete sub-results that are obsolete now.
+		delete leftStates;
+		delete rightStates;
+
+		LOG4CPLUS_INFO(logger, "Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states.");
+		LOG4CPLUS_INFO(logger, "Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states.");
+		storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
+		LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+
+		// Create resulting vector.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+
+		// Only try to solve system if there are states for which the probability is unknown.
+		uint_fast64_t maybeStatesSetBitCount = maybeStates.getNumberOfSetBits();
+		if (maybeStatesSetBitCount > 0) {
+			// First, we can eliminate the rows and columns from the original transition probability matrix for states
+			// whose probabilities are already known.
+			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates, *this->getModel().getNondeterministicChoiceIndices());
+
+			// Get the "new" nondeterministic choice indices for the submatrix.
+			std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
+
+			// Create vector for results for maybe states.
+			std::vector<Type> x(maybeStatesSetBitCount);
+
+			// Prepare the right-hand side of the equation system. For entry i this corresponds to
+			// the accumulated probability of going from state i to some 'yes' state.
+			std::vector<Type> b(submatrix->getRowCount());
+			this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, *this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, &b);
+
+			// Solve the corresponding system of equations.
+			this->solveEquationSystem(*submatrix, x, b, *subNondeterministicChoiceIndices);
+
+			// Set values of resulting vector according to result.
+			storm::utility::setVectorValues<Type>(result, maybeStates, x);
+
+			// Delete temporary matrix.
+			delete submatrix;
+		}
+
+		// Set values of resulting vector that are known exactly.
+		storm::utility::setVectorValues<Type>(result, statesWithProbability0, storm::utility::constGetZero<Type>());
+		storm::utility::setVectorValues<Type>(result, statesWithProbability1, storm::utility::constGetOne<Type>());
+
+		return result;
+	}
 
 	/*!
 	 * The check method for a path formula with an Instantaneous Reward operator node as root in its
@@ -255,7 +351,21 @@ public:
 	 * @param formula The Instantaneous Reward formula to check
 	 * @returns for each state the reward that the instantaneous reward yields
 	 */
-	virtual std::vector<Type>* checkInstantaneousReward(const storm::formula::InstantaneousReward<Type>& formula, bool qualitative) const = 0;
+	virtual std::vector<Type>* checkInstantaneousReward(const storm::formula::InstantaneousReward<Type>& formula, bool qualitative) const {
+		// Only compute the result if the model has a state-based reward model.
+		if (!this->getModel().hasStateRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing (state-based) reward model for formula.");
+			throw storm::exceptions::InvalidPropertyException() << "Missing (state-based) reward model for formula.";
+		}
+
+		// Initialize result to state rewards of the model.
+		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
+
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, nullptr, formula.getBound());
+
+		// Return result.
+		return result;
+	}
 
 	/*!
 	 * The check method for a path formula with a Cumulative Reward operator node as root in its
@@ -264,7 +374,32 @@ public:
 	 * @param formula The Cumulative Reward formula to check
 	 * @returns for each state the reward that the cumulative reward yields
 	 */
-	virtual std::vector<Type>* checkCumulativeReward(const storm::formula::CumulativeReward<Type>& formula, bool qualitative) const = 0;
+	virtual std::vector<Type>* checkCumulativeReward(const storm::formula::CumulativeReward<Type>& formula, bool qualitative) const {
+		// Only compute the result if the model has at least one reward model.
+		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing reward model for formula.");
+			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
+		}
+
+		// Compute the reward vector to add in each step based on the available reward models.
+		std::vector<Type>* totalRewardVector = nullptr;
+		if (this->getModel().hasTransitionRewards()) {
+			totalRewardVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
+			if (this->getModel().hasStateRewards()) {
+				gmm::add(*this->getModel().getStateRewardVector(), *totalRewardVector);
+			}
+		} else {
+			totalRewardVector = new std::vector<Type>(*this->getModel().getStateRewardVector());
+		}
+
+		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
+
+		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, totalRewardVector, formula.getBound());
+
+		// Delete temporary variables and return result.
+		delete totalRewardVector;
+		return result;
+	}
 
 	/*!
 	 * The check method for a path formula with a Reachability Reward operator node as root in its
@@ -273,10 +408,118 @@ public:
 	 * @param formula The Reachbility Reward formula to check
 	 * @returns for each state the reward that the reachability reward yields
 	 */
-	virtual std::vector<Type>* checkReachabilityReward(const storm::formula::ReachabilityReward<Type>& formula, bool qualitative) const = 0;
+	virtual std::vector<Type>* checkReachabilityReward(const storm::formula::ReachabilityReward<Type>& formula, bool qualitative) const {
+		// Only compute the result if the model has at least one reward model.
+		if (!this->getModel().hasStateRewards() && !this->getModel().hasTransitionRewards()) {
+			LOG4CPLUS_ERROR(logger, "Missing reward model for formula. Skipping formula");
+			throw storm::exceptions::InvalidPropertyException() << "Missing reward model for formula.";
+		}
+
+		// Determine the states for which the target predicate holds.
+		storm::storage::BitVector* targetStates = formula.getChild().check(*this);
+
+		// Determine which states have a reward of infinity by definition.
+		storm::storage::BitVector infinityStates(this->getModel().getNumberOfStates());
+		storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true);
+		if (this->minimumOperatorStack.top()) {
+			storm::utility::GraphAnalyzer::performProb1A(this->getModel(), trueStates, *targetStates, &infinityStates);
+		} else {
+			storm::utility::GraphAnalyzer::performProb1E(this->getModel(), trueStates, *targetStates, &infinityStates);
+		}
+		infinityStates.complement();
+
+		LOG4CPLUS_INFO(logger, "Found " << infinityStates.getNumberOfSetBits() << " 'infinity' states.");
+		LOG4CPLUS_INFO(logger, "Found " << targetStates->getNumberOfSetBits() << " 'target' states.");
+		storm::storage::BitVector maybeStates = ~(*targetStates) & ~infinityStates;
+		LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
+
+		// Create resulting vector.
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+
+		// Check whether there are states for which we have to compute the result.
+		const int maybeStatesSetBitCount = maybeStates.getNumberOfSetBits();
+		if (maybeStatesSetBitCount > 0) {
+			// First, we can eliminate the rows and columns from the original transition probability matrix for states
+			// whose probabilities are already known.
+			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates, *this->getModel().getNondeterministicChoiceIndices());
+
+			// Get the "new" nondeterministic choice indices for the submatrix.
+			std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
+
+			// Create vector for results for maybe states.
+			std::vector<Type> x(maybeStatesSetBitCount);
+
+			// Prepare the right-hand side of the equation system. For entry i this corresponds to
+			// the accumulated probability of going from state i to some 'yes' state.
+			std::vector<Type> b(submatrix->getRowCount());
+
+			if (this->getModel().hasTransitionRewards()) {
+				// If a transition-based reward model is available, we initialize the right-hand
+				// side to the vector resulting from summing the rows of the pointwise product
+				// of the transition probability matrix and the transition reward matrix.
+				std::vector<Type>* pointwiseProductRowSumVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
+				storm::utility::selectVectorValues(&b, maybeStates, *this->getModel().getNondeterministicChoiceIndices(), *pointwiseProductRowSumVector);
+				delete pointwiseProductRowSumVector;
+
+				if (this->getModel().hasStateRewards()) {
+					// If a state-based reward model is also available, we need to add this vector
+					// as well. As the state reward vector contains entries not just for the states
+					// that we still consider (i.e. maybeStates), we need to extract these values
+					// first.
+					std::vector<Type>* subStateRewards = new std::vector<Type>(b.size());
+					storm::utility::selectVectorValuesRepeatedly(subStateRewards, maybeStates, *this->getModel().getNondeterministicChoiceIndices(), *this->getModel().getStateRewardVector());
+					gmm::add(*subStateRewards, b);
+					delete subStateRewards;
+				}
+			} else {
+				// If only a state-based reward model is  available, we take this vector as the
+				// right-hand side. As the state reward vector contains entries not just for the
+				// states that we still consider (i.e. maybeStates), we need to extract these values
+				// first.
+				storm::utility::selectVectorValuesRepeatedly(&b, maybeStates, *this->getModel().getNondeterministicChoiceIndices(), *this->getModel().getStateRewardVector());
+			}
+
+			// Solve the corresponding system of equations.
+			this->solveEquationSystem(*submatrix, x, b, *subNondeterministicChoiceIndices);
+
+			// Set values of resulting vector according to result.
+			storm::utility::setVectorValues<Type>(result, maybeStates, x);
+			delete submatrix;
+		}
+
+		// Set values of resulting vector that are known exactly.
+		storm::utility::setVectorValues(result, *targetStates, storm::utility::constGetZero<Type>());
+		storm::utility::setVectorValues(result, infinityStates, storm::utility::constGetInfinity<Type>());
+
+		// Delete temporary storages and return result.
+		delete targetStates;
+		return result;
+	}
 
 protected:
 	mutable std::stack<bool> minimumOperatorStack;
+
+private:
+	virtual void performMatrixVectorMultiplication(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type>* summand = nullptr, uint_fast64_t repetitions = 1) const = 0;
+
+	virtual void solveEquationSystem(storm::storage::SparseMatrix<Type> const& matrix, std::vector<Type>& vector, std::vector<Type> const& b, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices) const = 0;
+
+
+	std::shared_ptr<std::vector<uint_fast64_t>> computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector constraint) const {
+		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
+		std::shared_ptr<std::vector<uint_fast64_t>> subNondeterministicChoiceIndices(new std::vector<uint_fast64_t>(constraint.getNumberOfSetBits() + 1));
+		uint_fast64_t currentRowCount = 0;
+		uint_fast64_t currentIndexCount = 1;
+		(*subNondeterministicChoiceIndices)[0] = 0;
+		for (auto index : constraint) {
+			(*subNondeterministicChoiceIndices)[currentIndexCount] = currentRowCount + (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index];
+			currentRowCount += (*nondeterministicChoiceIndices)[index + 1] - (*nondeterministicChoiceIndices)[index];
+			++currentIndexCount;
+		}
+		(*subNondeterministicChoiceIndices)[constraint.getNumberOfSetBits()] = currentRowCount;
+
+		return subNondeterministicChoiceIndices;
+	}
 };
 
 } //namespace modelChecker

test/functional/GmmxxMdpPrctModelCheckerTest.cpp

@@ -90,7 +90,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 
 	result = rewardFormula->check(mc);
 
-	ASSERT_LT(std::abs((*result)[0] - 7.3333272933959960938), s->get<double>("precision"));
+	ASSERT_LT(std::abs((*result)[0] - 7.3333294987678527832), s->get<double>("precision"));
 
 	delete rewardFormula;
 	delete result;
@@ -101,7 +101,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 
 	result = rewardFormula->check(mc);
 
-	ASSERT_LT(std::abs((*result)[0] - 7.3333272933959960938), s->get<double>("precision"));
+	ASSERT_LT(std::abs((*result)[0] - 7.3333294987678527832), s->get<double>("precision"));
 
 	delete rewardFormula;
 	delete result;
@@ -120,7 +120,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 
 	result = rewardFormula->check(stateRewardModelChecker);
 
-	ASSERT_LT(std::abs((*result)[0] - 7.3333272933959960938), s->get<double>("precision"));
+	ASSERT_LT(std::abs((*result)[0] - 7.3333294987678527832), s->get<double>("precision"));
 
 	delete rewardFormula;
 	delete result;
@@ -131,7 +131,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 
 	result = rewardFormula->check(stateRewardModelChecker);
 
-	ASSERT_LT(std::abs((*result)[0] - 7.3333272933959960938), s->get<double>("precision"));
+	ASSERT_LT(std::abs((*result)[0] - 7.3333294987678527832), s->get<double>("precision"));
 
 	delete rewardFormula;
 	delete result;
@@ -150,7 +150,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 
 	result = rewardFormula->check(stateAndTransitionRewardModelChecker);
 
-	ASSERT_LT(std::abs((*result)[0] - (2 * 7.3333272933959960938)), s->get<double>("precision"));
+	ASSERT_LT(std::abs((*result)[0] - (2 * 7.3333294987678527832)), s->get<double>("precision"));
 
 	delete rewardFormula;
 	delete result;
@@ -161,7 +161,7 @@ TEST(GmmxxMdpPrctModelCheckerTest, Dice) {
 
 	result = rewardFormula->check(stateAndTransitionRewardModelChecker);
 
-	ASSERT_LT(std::abs((*result)[0] - (2 * 7.3333272933959960938)), s->get<double>("precision"));
+	ASSERT_LT(std::abs((*result)[0] - (2 * 7.3333294987678527832)), s->get<double>("precision"));
 
 	delete rewardFormula;
 	delete result;