diff --git a/examples/mdp/two_dice/two_dice.nm b/examples/mdp/two_dice/two_dice.nm
new file mode 100644
index 000000000..e1bf34aea
--- /dev/null
+++ b/examples/mdp/two_dice/two_dice.nm
@@ -0,0 +1,40 @@
+// sum of two dice as the asynchronous parallel composition of
+// two copies of Knuth's model of a fair die using only fair coins
+
+mdp
+
+module die1
+
+	// local state
+	s1 : [0..7] init 0;
+	// value of the dice
+	d1 : [0..6] init 0;
+
+	[] s1=0 -> 0.5 : (s1'=1) + 0.5 : (s1'=2);
+	[] s1=1 -> 0.5 : (s1'=3) + 0.5 : (s1'=4);
+	[] s1=2 -> 0.5 : (s1'=5) + 0.5 : (s1'=6);
+	[] s1=3 -> 0.5 : (s1'=1) + 0.5 : (s1'=7) & (d1'=1);
+	[] s1=4 -> 0.5 : (s1'=7) & (d1'=2) + 0.5 : (s1'=7) & (d1'=3);
+	[] s1=5 -> 0.5 : (s1'=7) & (d1'=4) + 0.5 : (s1'=7) & (d1'=5);
+	[] s1=6 -> 0.5 : (s1'=2) + 0.5 : (s1'=7) & (d1'=6);
+	[] s1=7 & s2=7 -> (s1'=7);
+endmodule
+
+module die2 = die1 [ s1=s2, s2=s1, d1=d2 ] endmodule
+
+rewards "coinflips"
+	[] s1<7 | s2<7 : 1;
+endrewards
+
+label "done" = s1=7 & s2=7;
+label "two" = s1=7 & s2=7 & d1+d2=2;
+label "three" = s1=7 & s2=7 & d1+d2=3;
+label "four" = s1=7 & s2=7 & d1+d2=4;
+label "five" = s1=7 & s2=7 & d1+d2=5;
+label "six" = s1=7 & s2=7 & d1+d2=6;
+label "seven" = s1=7 & s2=7 & d1+d2=7;
+label "eight" = s1=7 & s2=7 & d1+d2=8;
+label "nine" = s1=7 & s2=7 & d1+d2=9;
+label "ten" = s1=7 & s2=7 & d1+d2=10;
+label "eleven" = s1=7 & s2=7 & d1+d2=11;
+label "twelve" = s1=7 & s2=7 & d1+d2=12;
diff --git a/src/modelchecker/GmmxxDtmcPrctlModelChecker.h b/src/modelchecker/GmmxxDtmcPrctlModelChecker.h
index a5a3b567d..9213df3fa 100644
--- a/src/modelchecker/GmmxxDtmcPrctlModelChecker.h
+++ b/src/modelchecker/GmmxxDtmcPrctlModelChecker.h
@@ -225,7 +225,7 @@ public:
 			totalRewardVector = new std::vector<Type>(*this->getModel().getStateRewardVector());
 		}
 
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
 
 		// Now perform matrix-vector multiplication as long as we meet the bound of the formula.
 		std::vector<Type>* swap = nullptr;
@@ -300,7 +300,7 @@ public:
 					// 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::setVectorValues(subStateRewards, maybeStates, *this->getModel().getStateRewardVector());
+					storm::utility::selectVectorValues(subStateRewards, maybeStates, *this->getModel().getStateRewardVector());
 					gmm::add(*subStateRewards, *b);
 					delete subStateRewards;
 				}
@@ -309,7 +309,7 @@ public:
 				// 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::setVectorValues(b, maybeStates, *this->getModel().getStateRewardVector());
+				storm::utility::selectVectorValues(b, maybeStates, *this->getModel().getStateRewardVector());
 			}
 
 			// Solve the corresponding system of linear equations.
diff --git a/src/modelchecker/GmmxxMdpPrctlModelChecker.h b/src/modelchecker/GmmxxMdpPrctlModelChecker.h
index e22da7263..a329fc303 100644
--- a/src/modelchecker/GmmxxMdpPrctlModelChecker.h
+++ b/src/modelchecker/GmmxxMdpPrctlModelChecker.h
@@ -197,24 +197,32 @@ public:
 			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.
-		std::vector<Type>* swap = nullptr;
-		std::vector<Type>* tmpResult = new std::vector<Type>(this->getModel().getNumberOfStates());
 		for (uint_fast64_t i = 0; i < formula.getBound(); ++i) {
-			gmm::mult(*gmmxxMatrix, *result, *tmpResult);
-			swap = tmpResult;
-			tmpResult = result;
-			result = swap;
+			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 tmpResult;
 		delete gmmxxMatrix;
 		return result;
 	}
@@ -226,6 +234,10 @@ public:
 			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());
 
@@ -240,23 +252,27 @@ public:
 			totalRewardVector = new std::vector<Type>(*this->getModel().getStateRewardVector());
 		}
 
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+		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.
-		std::vector<Type>* swap = nullptr;
-		std::vector<Type>* tmpResult = new std::vector<Type>(this->getModel().getNumberOfStates());
 		for (uint_fast64_t i = 0; i < formula.getBound(); ++i) {
-			gmm::mult(*gmmxxMatrix, *result, *tmpResult);
-			swap = tmpResult;
-			tmpResult = result;
-			result = swap;
+			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 tmpResult;
 		delete gmmxxMatrix;
 		delete totalRewardVector;
 		return result;
@@ -275,39 +291,52 @@ public:
 		// 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);
-		// TODO: just commented out to make it compile
-		//storm::utility::GraphAnalyzer::performProb1(this->getModel(), trueStates, *targetStates, &infinityStates);
+		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.
-		storm::storage::BitVector maybeStates = ~(*targetStates) & ~infinityStates;
 		const int maybeStatesSetBitCount = maybeStates.getNumberOfSetBits();
 		if (maybeStatesSetBitCount > 0) {
-			// Now we can eliminate the rows and columns from the original transition probability matrix.
-			storm::storage::SparseMatrix<Type>* submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates);
-			// Converting the matrix from the fixpoint notation to the form needed for the equation
-			// system. That is, we go from x = A*x + b to (I-A)x = b.
-			submatrix->convertToEquationSystem();
+			// 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 solvers.
+			// Transform the submatrix to the gmm++ format to use its capabilities.
 			gmm::csr_matrix<Type>* gmmxxMatrix = storm::adapters::GmmxxAdapter::toGmmxxSparseMatrix<Type>(*submatrix);
-			delete submatrix;
 
-			// Initialize the x vector with 1 for each element. This is the initial guess for
-			// the iterative solvers.
-			std::vector<Type> x(maybeStatesSetBitCount, storm::utility::constGetOne<Type>());
+			// 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;
 
-			// Prepare the right-hand side of the equation system.
-			std::vector<Type>* b = new std::vector<Type>(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, *nondeterministicChoiceIndices, *pointwiseProductRowSumVector);
 				delete pointwiseProductRowSumVector;
 
 				if (this->getModel().hasStateRewards()) {
@@ -315,9 +344,9 @@ 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::setVectorValues(subStateRewards, maybeStates, *this->getModel().getStateRewardVector());
-					gmm::add(*subStateRewards, *b);
+					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 {
@@ -325,19 +354,16 @@ public:
 				// 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::setVectorValues(b, maybeStates, *this->getModel().getStateRewardVector());
+				storm::utility::selectVectorValuesRepeatedly(&b, maybeStates, *nondeterministicChoiceIndices, *this->getModel().getStateRewardVector());
 			}
 
-			// Solve the corresponding system of linear equations.
-			// TODO: just commented out to make it compile
-			// this->solveEquationSystem(*gmmxxMatrix, x, *b);
+			// 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 right-hand side.
+			storm::utility::setVectorValues<Type>(result, maybeStates, *x);
+			delete x;
 			delete gmmxxMatrix;
-			delete b;
 		}
 
 		// Set values of resulting vector that are known exactly.
@@ -400,6 +426,11 @@ private:
 			++iterations;
 		}
 
+		if (iterations % 2 == 1) {
+			delete x;
+		} else {
+			delete newX;
+		}
 		delete temporaryResult;
 
 		// Check if the solver converged and issue a warning otherwise.
diff --git a/src/models/Dtmc.h b/src/models/Dtmc.h
index eefcc2ab5..7f0a0bb89 100644
--- a/src/models/Dtmc.h
+++ b/src/models/Dtmc.h
@@ -50,10 +50,11 @@ public:
 			throw storm::exceptions::InvalidArgumentException() << "Probability matrix is invalid.";
 		}
 		if (this->getTransitionRewardMatrix() != nullptr) {
-			if (!this->getTransitionRewardMatrix()->isSubmatrixOf(*(this->getTransitionMatrix()))) {
-				LOG4CPLUS_ERROR(logger, "Transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
-				throw storm::exceptions::InvalidArgumentException() << "There are transition rewards for nonexistent transitions.";
-			}
+			// FIXME: This is temporarily commented out, because the checking routine does not work properly.
+			//if (!this->getTransitionRewardMatrix()->isSubmatrixOf(*(this->getTransitionMatrix()))) {
+			//	LOG4CPLUS_ERROR(logger, "Transition reward matrix is not a submatrix of the transition matrix, i.e. there are rewards for transitions that do not exist.");
+			//	throw storm::exceptions::InvalidArgumentException() << "There are transition rewards for nonexistent transitions.";
+			//}
 		}
 	}
 
diff --git a/src/parser/DeterministicModelParser.cpp b/src/parser/DeterministicModelParser.cpp
index 2372cde3e..2300ad74d 100644
--- a/src/parser/DeterministicModelParser.cpp
+++ b/src/parser/DeterministicModelParser.cpp
@@ -45,7 +45,7 @@ DeterministicModelParser::DeterministicModelParser(std::string const & transitio
 		this->stateRewards = srp.getStateRewards();
 	}
 	if (transitionRewardFile != "") {
-		storm::parser::DeterministicSparseTransitionParser trp(transitionRewardFile);
+		storm::parser::DeterministicSparseTransitionParser trp(transitionRewardFile, false, true);
 		this->transitionRewards = trp.getMatrix();
 	}
 	this->dtmc = nullptr;
diff --git a/src/parser/DeterministicSparseTransitionParser.cpp b/src/parser/DeterministicSparseTransitionParser.cpp
index a91c61525..031087b88 100644
--- a/src/parser/DeterministicSparseTransitionParser.cpp
+++ b/src/parser/DeterministicSparseTransitionParser.cpp
@@ -44,12 +44,14 @@ namespace parser {
  *	@param buf Data to scan. Is expected to be some char array.
  *	@param maxnode Is set to highest id of all nodes.
  */
-uint_fast64_t DeterministicSparseTransitionParser::firstPass(char* buf, int_fast64_t& maxnode) {
+uint_fast64_t DeterministicSparseTransitionParser::firstPass(char* buf, int_fast64_t& maxnode, bool rewardFile) {
 	uint_fast64_t nonZeroEntryCount = 0;
 	/*
 	 *	Check file header and extract number of transitions.
 	 */
-	buf = strchr(buf, '\n') + 1;  // skip format hint
+	if (!rewardFile) {
+		buf = strchr(buf, '\n') + 1;  // skip format hint
+	}
 
     /*
      * Check all transitions for non-zero diagonal entries and deadlock states.
@@ -66,24 +68,26 @@ uint_fast64_t DeterministicSparseTransitionParser::firstPass(char* buf, int_fast
             row = checked_strtol(buf, &buf);
             col = checked_strtol(buf, &buf);
 
-			if (lastRow != row) {
-				if ((lastRow != -1) && (!rowHadDiagonalEntry)) {
-					++nonZeroEntryCount;
-					rowHadDiagonalEntry = true;
+            if (!rewardFile) {
+				if (lastRow != row) {
+					if ((lastRow != -1) && (!rowHadDiagonalEntry)) {
+						++nonZeroEntryCount;
+						rowHadDiagonalEntry = true;
+					}
+					for (int_fast64_t skippedRow = lastRow + 1; skippedRow < row; ++skippedRow) {
+						++nonZeroEntryCount;
+					}
+					lastRow = row;
+					rowHadDiagonalEntry = false;
 				}
-				for (int_fast64_t skippedRow = lastRow + 1; skippedRow < row; ++skippedRow) {
+
+				if (col == row) {
+					rowHadDiagonalEntry = true;
+				} else if (col > row && !rowHadDiagonalEntry) {
+					rowHadDiagonalEntry = true;
 					++nonZeroEntryCount;
 				}
-				lastRow = row;
-				rowHadDiagonalEntry = false;
-			}
-
-			if (col == row) {
-				rowHadDiagonalEntry = true;
-			} else if (col > row && !rowHadDiagonalEntry) {
-				rowHadDiagonalEntry = true;
-				++nonZeroEntryCount;
-			}
+            }
 
             /*
              *      Check if one is larger than the current maximum id.
@@ -105,7 +109,7 @@ uint_fast64_t DeterministicSparseTransitionParser::firstPass(char* buf, int_fast
             buf = trimWhitespaces(buf);
     }
 
-    if (!rowHadDiagonalEntry) {
+    if (!rowHadDiagonalEntry && !rewardFile) {
     	++nonZeroEntryCount;
     }
 
@@ -122,7 +126,7 @@ uint_fast64_t DeterministicSparseTransitionParser::firstPass(char* buf, int_fast
  *	@return a pointer to the created sparse matrix.
  */
 
-DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::string const &filename, bool insertDiagonalEntriesIfMissing)
+DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::string const &filename, bool insertDiagonalEntriesIfMissing, bool rewardFile)
 	: matrix(nullptr) {
 	/*
 	 *	Enforce locale where decimal point is '.'.
@@ -139,7 +143,7 @@ DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::st
 	 *	Perform first pass, i.e. count entries that are not zero.
 	 */
 	int_fast64_t maxStateId;
-	uint_fast64_t nonZeroEntryCount = this->firstPass(file.data, maxStateId);
+	uint_fast64_t nonZeroEntryCount = this->firstPass(file.data, maxStateId, rewardFile);
 
 	LOG4CPLUS_INFO(logger, "First pass on " << filename << " shows " << nonZeroEntryCount << " NonZeros.");
 
@@ -160,7 +164,9 @@ DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::st
 	/*
 	 *	Read file header, extract number of states.
 	 */
-	buf = strchr(buf, '\n') + 1;  // skip format hint
+	if (!rewardFile) {
+		buf = strchr(buf, '\n') + 1;  // skip format hint
+	}
 
 	/*
 	 *	Creating matrix here.
@@ -198,7 +204,7 @@ DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::st
 				if (insertDiagonalEntriesIfMissing) {
 					this->matrix->addNextValue(lastRow, lastRow, storm::utility::constGetZero<double>());
 					LOG4CPLUS_DEBUG(logger, "While parsing " << filename << ": state " << lastRow << " has no transition to itself. Inserted a 0-transition. (1)");
-				} else {
+				} else if (!rewardFile) {
 					LOG4CPLUS_WARN(logger, "Warning while parsing " << filename << ": state " << lastRow << " has no transition to itself.");
 				}
 				// No increment for lastRow
@@ -206,11 +212,11 @@ DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::st
 			}
 			for (int_fast64_t skippedRow = lastRow + 1; skippedRow < row; ++skippedRow) {
 				hadDeadlocks = true;
-				if (fixDeadlocks) {
+				if (fixDeadlocks && !rewardFile) {
 					this->matrix->addNextValue(skippedRow, skippedRow, storm::utility::constGetOne<double>());
 					rowHadDiagonalEntry = true;
 					LOG4CPLUS_WARN(logger, "Warning while parsing " << filename << ": state " << skippedRow << " has no outgoing transitions. A self-loop was inserted.");
-				} else {
+				} else if (!rewardFile) {
 					LOG4CPLUS_ERROR(logger, "Error while parsing " << filename << ": state " << skippedRow << " has no outgoing transitions.");
 					// FIXME Why no exception at this point? This will break the App.
 				}
@@ -223,10 +229,10 @@ DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::st
 			rowHadDiagonalEntry = true;
 		} else if (col > row && !rowHadDiagonalEntry) {
 			rowHadDiagonalEntry = true;
-			if (insertDiagonalEntriesIfMissing) {
+			if (insertDiagonalEntriesIfMissing && !rewardFile) {
 				this->matrix->addNextValue(row, row, storm::utility::constGetZero<double>());
 				LOG4CPLUS_DEBUG(logger, "While parsing " << filename << ": state " << row << " has no transition to itself. Inserted a 0-transition. (2)");
-			} else {
+			} else if (!rewardFile) {
 				LOG4CPLUS_WARN(logger, "Warning while parsing " << filename << ": state " << row << " has no transition to itself.");
 			}
 		}
@@ -236,10 +242,10 @@ DeterministicSparseTransitionParser::DeterministicSparseTransitionParser(std::st
 	}
 
 	if (!rowHadDiagonalEntry) {
-		if (insertDiagonalEntriesIfMissing) {
+		if (insertDiagonalEntriesIfMissing && !rewardFile) {
 			this->matrix->addNextValue(lastRow, lastRow, storm::utility::constGetZero<double>());
 			LOG4CPLUS_DEBUG(logger, "While parsing " << filename << ": state " << lastRow << " has no transition to itself. Inserted a 0-transition. (3)");
-		} else {
+		} else if (!rewardFile) {
 			LOG4CPLUS_WARN(logger, "Warning while parsing " << filename << ": state " << lastRow << " has no transition to itself.");
 		}
 	}
diff --git a/src/parser/DeterministicSparseTransitionParser.h b/src/parser/DeterministicSparseTransitionParser.h
index 2407f3398..a023debb7 100644
--- a/src/parser/DeterministicSparseTransitionParser.h
+++ b/src/parser/DeterministicSparseTransitionParser.h
@@ -17,7 +17,7 @@ namespace parser {
  */
 class DeterministicSparseTransitionParser : public Parser {
 	public:
-		DeterministicSparseTransitionParser(std::string const &filename, bool insertDiagonalEntriesIfMissing = true);
+		DeterministicSparseTransitionParser(std::string const &filename, bool insertDiagonalEntriesIfMissing = true, bool rewardFile = false);
 		
 		std::shared_ptr<storm::storage::SparseMatrix<double>> getMatrix() {
 			return this->matrix;
@@ -26,7 +26,7 @@ class DeterministicSparseTransitionParser : public Parser {
 	private:
 		std::shared_ptr<storm::storage::SparseMatrix<double>> matrix;
 		
-		uint_fast64_t firstPass(char* buf, int_fast64_t &maxnode);
+		uint_fast64_t firstPass(char* buf, int_fast64_t &maxnode, bool rewardFile);
 	
 };
 		
diff --git a/src/parser/NondeterministicModelParser.cpp b/src/parser/NondeterministicModelParser.cpp
index 87eb3cc57..303789a99 100644
--- a/src/parser/NondeterministicModelParser.cpp
+++ b/src/parser/NondeterministicModelParser.cpp
@@ -30,7 +30,7 @@ namespace parser {
 NondeterministicModelParser::NondeterministicModelParser(std::string const & transitionSystemFile, std::string const & labelingFile,
 		std::string const & stateRewardFile, std::string const & transitionRewardFile) {
 	storm::parser::NondeterministicSparseTransitionParser tp(transitionSystemFile);
-	uint_fast64_t stateCount = tp.getMatrix()->getRowCount();
+	uint_fast64_t stateCount = tp.getMatrix()->getColumnCount();
 
 	storm::parser::AtomicPropositionLabelingParser lp(stateCount, labelingFile);
 	if (stateRewardFile != "") {
@@ -38,7 +38,7 @@ NondeterministicModelParser::NondeterministicModelParser(std::string const & tra
 		this->stateRewards = srp.getStateRewards();
 	}
 	if (transitionRewardFile != "") {
-		storm::parser::NondeterministicSparseTransitionParser trp(transitionRewardFile);
+		storm::parser::NondeterministicSparseTransitionParser trp(transitionRewardFile, true, tp.getRowMapping());
 		this->transitionRewardMatrix = trp.getMatrix();
 	}
 
diff --git a/src/parser/NondeterministicSparseTransitionParser.cpp b/src/parser/NondeterministicSparseTransitionParser.cpp
index 66089c86f..1019240c7 100644
--- a/src/parser/NondeterministicSparseTransitionParser.cpp
+++ b/src/parser/NondeterministicSparseTransitionParser.cpp
@@ -49,11 +49,13 @@ namespace parser {
  *	@param maxnode Is set to highest id of all nodes.
  *	@return The number of non-zero elements.
  */
-uint_fast64_t NondeterministicSparseTransitionParser::firstPass(char* buf, uint_fast64_t& choices, int_fast64_t& maxnode) {
+uint_fast64_t NondeterministicSparseTransitionParser::firstPass(char* buf, uint_fast64_t& choices, int_fast64_t& maxnode, bool rewardFile, std::shared_ptr<std::vector<uint_fast64_t>> const nondeterministicChoiceIndices) {
 	/*
 	 *	Check file header and extract number of transitions.
 	 */
-	buf = strchr(buf, '\n') + 1;  // skip format hint
+	if (!rewardFile) {
+		buf = strchr(buf, '\n') + 1;  // skip format hint
+	}
 
 	/*
 	 *	Read all transitions.
@@ -78,15 +80,37 @@ uint_fast64_t NondeterministicSparseTransitionParser::firstPass(char* buf, uint_
 			maxnode = source;
 		}
 
-		// If we have skipped some states, we need to reserve the space for the self-loop insertion
-		// in the second pass.
-		if (source > lastsource + 1) {
-			nonzero += source - lastsource - 1;
-			choices += source - lastsource - 1;
-		} else if (source != lastsource || choice != lastchoice) {
-			// If we have switched the source state or the nondeterministic choice, we need to
-			// reserve one row more.
-			++choices;
+		if (rewardFile) {
+			// If we have switched the source state, we possibly need to insert the rows of the last
+			// last source state.
+			if (source != lastsource && lastsource != -1) {
+				choices += lastchoice - ((*nondeterministicChoiceIndices)[lastsource + 1] - (*nondeterministicChoiceIndices)[lastsource] - 1);
+			}
+
+			// If we skipped some states, we need to reserve empty rows for all their nondeterministic
+			// choices.
+			for (uint_fast64_t i = lastsource + 1; i < source; ++i) {
+				choices += ((*nondeterministicChoiceIndices)[i + 1] - (*nondeterministicChoiceIndices)[i]);
+			}
+
+			// If we advanced to the next state, but skipped some choices, we have to reserve rows
+			// for them
+			if (source != lastsource) {
+				choices += choice + 1;
+			} else if (choice != lastchoice) {
+				choices += choice - lastchoice;
+			}
+		} else {
+			// If we have skipped some states, we need to reserve the space for the self-loop insertion
+			// in the second pass.
+			if (source > lastsource + 1) {
+				nonzero += source - lastsource - 1;
+				choices += source - lastsource - 1;
+			} else if (source != lastsource || choice != lastchoice) {
+				// If we have switched the source state or the nondeterministic choice, we need to
+				// reserve one row more.
+				++choices;
+			}
 		}
 
 		// Read target and check if we encountered a state index that is bigger than all previously
@@ -127,6 +151,17 @@ uint_fast64_t NondeterministicSparseTransitionParser::firstPass(char* buf, uint_
 		buf = trimWhitespaces(buf);
 	}
 
+	if (rewardFile) {
+		// If not all rows were filled for the last state, we need to insert them.
+		choices += lastchoice - ((*nondeterministicChoiceIndices)[lastsource + 1] - (*nondeterministicChoiceIndices)[lastsource] - 1);
+
+		// If we skipped some states, we need to reserve empty rows for all their nondeterministic
+		// choices.
+		for (uint_fast64_t i = lastsource + 1; i < nondeterministicChoiceIndices->size() - 1; ++i) {
+			choices += ((*nondeterministicChoiceIndices)[i + 1] - (*nondeterministicChoiceIndices)[i]);
+		}
+	}
+
 	return nonzero;
 }
 
@@ -140,7 +175,7 @@ uint_fast64_t NondeterministicSparseTransitionParser::firstPass(char* buf, uint_
  *	@return a pointer to the created sparse matrix.
  */
 
-NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(std::string const &filename)
+NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(std::string const &filename, bool rewardFile, std::shared_ptr<std::vector<uint_fast64_t>> const nondeterministicChoiceIndices)
 	: matrix(nullptr) {
 	/*
 	 *	Enforce locale where decimal point is '.'.
@@ -158,7 +193,7 @@ NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(s
 	 */
 	int_fast64_t maxnode;
 	uint_fast64_t choices;
-	uint_fast64_t nonzero = this->firstPass(file.data, choices, maxnode);
+	uint_fast64_t nonzero = this->firstPass(file.data, choices, maxnode, rewardFile, nondeterministicChoiceIndices);
 
 	/*
 	 *	If first pass returned zero, the file format was wrong.
@@ -175,9 +210,11 @@ NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(s
 	 */
 
 	/*
-	 *	Read file header, ignore values within.
+	 *	Skip file header.
 	 */
-	buf = strchr(buf, '\n') + 1;  // skip format hint
+	if (!rewardFile) {
+		buf = strchr(buf, '\n') + 1;  // skip format hint
+	}
 
 	/*
 	 *	Create and initialize matrix.
@@ -216,33 +253,54 @@ NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(s
 		source = checked_strtol(buf, &buf);
 		choice = checked_strtol(buf, &buf);
 
-		// Increase line count if we have either finished reading the transitions of a certain state
-		// or we have finished reading one nondeterministic choice of a state.
-		if ((source != lastsource || choice != lastchoice)) {
-			++curRow;
-		}
+		if (rewardFile) {
+			// If we have switched the source state, we possibly need to insert the rows of the last
+			// last source state.
+			if (source != lastsource && lastsource != -1) {
+				curRow += lastchoice - ((*nondeterministicChoiceIndices)[lastsource + 1] - (*nondeterministicChoiceIndices)[lastsource] - 1);
+			}
 
-		/*
-		 *	Check if we have skipped any source node, i.e. if any node has no
-		 *	outgoing transitions. If so, insert a self-loop.
-		 *	Also add self-loops to rowMapping.
-		 */
-		for (int_fast64_t node = lastsource + 1; node < source; node++) {
-			hadDeadlocks = true;
-			if (fixDeadlocks) {
-				this->rowMapping->at(node) = curRow;
-				this->matrix->addNextValue(curRow, node, 1);
+			// If we skipped some states, we need to reserve empty rows for all their nondeterministic
+			// choices.
+			for (uint_fast64_t i = lastsource + 1; i < source; ++i) {
+				curRow += ((*nondeterministicChoiceIndices)[i + 1] - (*nondeterministicChoiceIndices)[i]);
+			}
+
+			// If we advanced to the next state, but skipped some choices, we have to reserve rows
+			// for them
+			if (source != lastsource) {
+				curRow += choice + 1;
+			} else if (choice != lastchoice) {
+				curRow += choice - lastchoice;
+			}
+		} else {
+			// Increase line count if we have either finished reading the transitions of a certain state
+			// or we have finished reading one nondeterministic choice of a state.
+			if ((source != lastsource || choice != lastchoice)) {
 				++curRow;
-				LOG4CPLUS_WARN(logger, "Warning while parsing " << filename << ": node " << node << " has no outgoing transitions. A self-loop was inserted.");
-			} else {
-				LOG4CPLUS_ERROR(logger, "Error while parsing " << filename << ": node " << node << " has no outgoing transitions.");
 			}
-		}
-		if (source != lastsource) {
 			/*
-			 *	Add this source to rowMapping, if this is the first choice we encounter for this state.
+			 *	Check if we have skipped any source node, i.e. if any node has no
+			 *	outgoing transitions. If so, insert a self-loop.
+			 *	Also add self-loops to rowMapping.
 			 */
-			this->rowMapping->at(source) = curRow;
+			for (int_fast64_t node = lastsource + 1; node < source; node++) {
+				hadDeadlocks = true;
+				if (!rewardFile && fixDeadlocks) {
+					this->rowMapping->at(node) = curRow;
+					this->matrix->addNextValue(curRow, node, 1);
+					++curRow;
+					LOG4CPLUS_WARN(logger, "Warning while parsing " << filename << ": node " << node << " has no outgoing transitions. A self-loop was inserted.");
+				} else if (!rewardFile) {
+					LOG4CPLUS_ERROR(logger, "Error while parsing " << filename << ": node " << node << " has no outgoing transitions.");
+				}
+			}
+			if (source != lastsource) {
+				/*
+				 *	Add this source to rowMapping, if this is the first choice we encounter for this state.
+				 */
+				this->rowMapping->at(source) = curRow;
+			}
 		}
 
 		// Read target and value and write it to the matrix.
@@ -265,7 +323,7 @@ NondeterministicSparseTransitionParser::NondeterministicSparseTransitionParser(s
 
 	this->rowMapping->at(maxnode+1) = curRow + 1;
 
-	if (!fixDeadlocks && hadDeadlocks) throw storm::exceptions::WrongFileFormatException() << "Some of the nodes had deadlocks. You can use --fix-deadlocks to insert self-loops on the fly.";
+	if (!fixDeadlocks && hadDeadlocks && !rewardFile) throw storm::exceptions::WrongFileFormatException() << "Some of the nodes had deadlocks. You can use --fix-deadlocks to insert self-loops on the fly.";
 
 	/*
 	 * Finalize matrix.
diff --git a/src/parser/NondeterministicSparseTransitionParser.h b/src/parser/NondeterministicSparseTransitionParser.h
index 209c0e7a0..79947dd93 100644
--- a/src/parser/NondeterministicSparseTransitionParser.h
+++ b/src/parser/NondeterministicSparseTransitionParser.h
@@ -19,7 +19,7 @@ namespace parser {
  */
 class NondeterministicSparseTransitionParser : public Parser {
 	public:
-		NondeterministicSparseTransitionParser(std::string const &filename);
+		NondeterministicSparseTransitionParser(std::string const &filename, bool rewardFile = false, std::shared_ptr<std::vector<uint_fast64_t>> const nondeterministicChoiceIndices = nullptr);
 		
 		inline std::shared_ptr<storm::storage::SparseMatrix<double>> getMatrix() const {
 			return this->matrix;
@@ -33,7 +33,7 @@ class NondeterministicSparseTransitionParser : public Parser {
 		std::shared_ptr<storm::storage::SparseMatrix<double>> matrix;
 		std::shared_ptr<std::vector<uint_fast64_t>> rowMapping;
 		
-		uint_fast64_t firstPass(char* buf, uint_fast64_t& choices, int_fast64_t& maxnode);
+		uint_fast64_t firstPass(char* buf, uint_fast64_t& choices, int_fast64_t& maxnode, bool rewardFile, std::shared_ptr<std::vector<uint_fast64_t>> const nondeterministicChoiceIndices);
 	
 };
 		
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index f45ec4f96..bde75c354 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -226,8 +226,8 @@ public:
 			throw storm::exceptions::InvalidStateException("Trying to finalize an uninitialized matrix.");
 		} else if (currentSize != nonZeroEntryCount) {
 			triggerErrorState();
-			LOG4CPLUS_ERROR(logger, "Trying to finalize a matrix that was initialized with more non-zero entries than given.");
-			throw storm::exceptions::InvalidStateException("Trying to finalize a matrix that was initialized with more non-zero entries than given.");
+			LOG4CPLUS_ERROR(logger, "Trying to finalize a matrix that was initialized with more non-zero entries than given (expected " << nonZeroEntryCount << " but got " << currentSize << " instead)");
+			throw storm::exceptions::InvalidStateException() << "Trying to finalize a matrix that was initialized with more non-zero entries than given (expected " << nonZeroEntryCount << " but got " << currentSize << " instead).";
 		} else {
 			// Fill in the missing entries in the row_indications array.
 			// (Can happen because of empty rows at the end.)
@@ -751,7 +751,7 @@ public:
 	 */
 	std::vector<T>* getPointwiseProductRowSumVector(storm::storage::SparseMatrix<T> const& otherMatrix) {
 		// Prepare result.
-		std::vector<T>* result = new std::vector<T>(rowCount);
+		std::vector<T>* result = new std::vector<T>(rowCount, storm::utility::constGetZero<T>());
 
 		// Iterate over all elements of the current matrix and either continue with the next element
 		// in case the given matrix does not have a non-zero element at this column position, or
@@ -764,7 +764,7 @@ public:
 					// If the precondition of this method (i.e. that the given matrix is a submatrix
 					// of the current one) was fulfilled, we know now that the two elements are in
 					// the same column, so we can multiply and add them to the row sum vector.
-					(*result)[row] += otherMatrix.valueStorage[element] * valueStorage[nextOtherElement];
+					(*result)[row] += otherMatrix.valueStorage[nextOtherElement] * valueStorage[element];
 					++nextOtherElement;
 				}
 			}
@@ -866,6 +866,7 @@ public:
 	 *	@return True iff this is a submatrix of matrix.
 	 */
 	bool isSubmatrixOf(SparseMatrix<T> const & matrix) const {
+		// FIXME: THIS DOES NOT IMPLEMENT WHAT IS PROMISED.
 		if (this->getRowCount() != matrix.getRowCount()) return false;
 		if (this->getColumnCount() != matrix.getColumnCount()) return false;
 
@@ -933,7 +934,7 @@ public:
 			result << i << "\t(\t";
 			uint_fast64_t currentRealIndex = 0;
 			while (currentRealIndex < colCount) {
-				if (currentRealIndex == columnIndications[nextIndex] && nextIndex < rowIndications[i + 1]) {
+				if (nextIndex < rowIndications[i + 1] && currentRealIndex == columnIndications[nextIndex]) {
 					result << valueStorage[nextIndex] << "\t";
 					++nextIndex;
 				} else {
diff --git a/src/storm.cpp b/src/storm.cpp
index d75982d34..f2b34cdbd 100644
--- a/src/storm.cpp
+++ b/src/storm.cpp
@@ -220,6 +220,14 @@ void testCheckingDice(storm::models::Mdp<double>& mdp) {
 	mc->check(*probFormula);
 	delete probFormula;
 
+	storm::formula::Ap<double>* doneFormula = new storm::formula::Ap<double>("done");
+	//storm::formula::ReachabilityReward<double>* reachabilityRewardFormula = new storm::formula::ReachabilityReward<double>(doneFormula);
+	storm::formula::InstantaneousReward<double>* reachabilityRewardFormula = new storm::formula::InstantaneousReward<double>(8);
+	storm::formula::RewardNoBoundOperator<double>* rewardFormula = new storm::formula::RewardNoBoundOperator<double>(reachabilityRewardFormula, true);
+
+	mc->check(*rewardFormula);
+	delete rewardFormula;
+
 	delete mc;
 }
 
@@ -271,8 +279,7 @@ void testChecking() {
 		// testCheckingDie(*dtmc);
 		// testCheckingCrowds(*dtmc);
 		// testCheckingSynchronousLeader(*dtmc, 4);
-	}
-	else if (parser.getType() == storm::models::MDP) {
+	} else if (parser.getType() == storm::models::MDP) {
 		std::shared_ptr<storm::models::Mdp<double>> mdp = parser.getModel<storm::models::Mdp<double>>();
 		mdp->printModelInformationToStream(std::cout);
 
diff --git a/src/utility/GraphAnalyzer.h b/src/utility/GraphAnalyzer.h
index c7325befd..f9b6b783a 100644
--- a/src/utility/GraphAnalyzer.h
+++ b/src/utility/GraphAnalyzer.h
@@ -262,6 +262,8 @@ public:
 			} else {
 				*currentStates = *nextStates;
 			}
+
+			delete nextStates;
 		}
 
 		*statesWithProbability1 = *currentStates;
@@ -407,6 +409,7 @@ public:
 			} else {
 				*currentStates = *nextStates;
 			}
+			delete nextStates;
 		}
 
 		*statesWithProbability1 = *currentStates;
diff --git a/src/utility/Vector.h b/src/utility/Vector.h
index 7e5224af4..fce8c18e1 100644
--- a/src/utility/Vector.h
+++ b/src/utility/Vector.h
@@ -51,6 +51,26 @@ void selectVectorValues(std::vector<T>* vector, const storm::storage::BitVector&
 	}
 }
 
+template<class T>
+void selectVectorValues(std::vector<T>* vector, const storm::storage::BitVector& positions, const std::vector<uint_fast64_t>& rowMapping, std::vector<T> const& values) {
+	uint_fast64_t oldPosition = 0;
+	for (auto position : positions) {
+		for (uint_fast64_t i = rowMapping[position]; i < rowMapping[position + 1]; ++i) {
+			(*vector)[oldPosition++] = values[i];
+		}
+	}
+}
+
+template<class T>
+void selectVectorValuesRepeatedly(std::vector<T>* vector, const storm::storage::BitVector& positions, const std::vector<uint_fast64_t>& rowMapping, std::vector<T> const& values) {
+	uint_fast64_t oldPosition = 0;
+	for (auto position : positions) {
+		for (uint_fast64_t i = rowMapping[position]; i < rowMapping[position + 1]; ++i) {
+			(*vector)[oldPosition++] = values[position];
+		}
+	}
+}
+
 template<class T>
 void subtractFromConstantOneVector(std::vector<T>* vector) {
 	for (auto it = vector->begin(); it != vector->end(); ++it) {
@@ -101,7 +121,7 @@ void reduceVectorMax(std::vector<T> const& source, std::vector<T>* target, std::
 template<class T>
 bool equalModuloPrecision(std::vector<T> const& vectorLeft, std::vector<T> const& vectorRight, T precision, bool relativeError) {
 	if (vectorLeft.size() != vectorRight.size()) {
-		LOG4CPLUS_ERROR(logger, "Length of vectors does not match and makes comparison impossible.");
+		LOG4CPLUS_ERROR(logger, "Lengths of vectors does not match and makes comparison impossible.");
 		throw storm::exceptions::InvalidArgumentException() << "Length of vectors does not match and makes comparison impossible.";
 	}