diff --git a/src/modelchecker/prctl/AbstractModelChecker.h b/src/modelchecker/prctl/AbstractModelChecker.h
index 2ac3fd241..b3ae46bd2 100644
--- a/src/modelchecker/prctl/AbstractModelChecker.h
+++ b/src/modelchecker/prctl/AbstractModelChecker.h
@@ -86,7 +86,8 @@ public:
 
 	/*!
 	 * Returns a pointer to the model checker object that is of the requested type as given by the template parameters.
-	 * @returns A pointer to the model checker object that is of the requested type as given by the template parameters.
+     *
+	 * @return A pointer to the model checker object that is of the requested type as given by the template parameters.
 	 * If the model checker is not of the requested type, type casting will fail and result in an exception.
 	 */
 	template <template <class T> class Target>
@@ -105,7 +106,7 @@ public:
 	 * Retrieves the model associated with this model checker as a constant reference to an object of the type given
 	 * by the template parameter.
 	 *
-	 * @returns A constant reference of the specified type to the model associated with this model checker. If the model
+	 * @return A constant reference of the specified type to the model associated with this model checker. If the model
 	 * is not of the requested type, type casting will fail and result in an exception.
 	 */
 	template <class Model>
@@ -118,6 +119,15 @@ public:
 			throw bc;
 		}
 	}
+        
+    /*!
+     * Retrieves the initial states of the model.
+     *
+     * @return A bit vector that represents the initial states of the model.
+     */
+    storm::storage::BitVector const& getInitialStates() const {
+        return model.getLabeledStates("init");
+    }
 
 	/*!
 	 * Checks the given abstract prctl formula on the model and prints the result (depending on the actual type of the formula)
@@ -148,7 +158,7 @@ public:
 			result = stateFormula.check(*this);
 			LOG4CPLUS_INFO(logger, "Result for initial states:");
 			std::cout << "Result for initial states:" << std::endl;
-			for (auto initialState : model.getLabeledStates("init")) {
+			for (auto initialState : this->getInitialStates()) {
 				LOG4CPLUS_INFO(logger, "\t" << initialState << ": " << (result->get(initialState) ? "satisfied" : "not satisfied"));
 				std::cout << "\t" << initialState << ": " << result->get(initialState) << std::endl;
 			}
@@ -178,7 +188,7 @@ public:
 			result = this->checkNoBoundOperator(noBoundFormula);
 			LOG4CPLUS_INFO(logger, "Result for initial states:");
 			std::cout << "Result for initial states:" << std::endl;
-			for (auto initialState : model.getLabeledStates("init")) {
+			for (auto initialState : this->getInitialStates()) {
 				LOG4CPLUS_INFO(logger, "\t" << initialState << ": " << (*result)[initialState]);
 				std::cout << "\t" << initialState << ": " << (*result)[initialState] << std::endl;
 			}
@@ -196,7 +206,7 @@ public:
 	 * Checks the given formula consisting of a single atomic proposition.
 	 *
 	 * @param formula The formula to be checked.
-	 * @returns The set of states satisfying the formula represented by a bit vector.
+	 * @return The set of states satisfying the formula represented by a bit vector.
 	 */
 	storm::storage::BitVector* checkAp(storm::property::prctl::Ap<Type> const& formula) const {
 		if (formula.getAp() == "true") {
@@ -217,7 +227,7 @@ public:
 	 * Checks the given formula that is a logical "and" of two formulae.
 	 *
 	 * @param formula The formula to be checked.
-	 * @returns The set of states satisfying the formula represented by a bit vector.
+	 * @return The set of states satisfying the formula represented by a bit vector.
 	 */
 	storm::storage::BitVector* checkAnd(storm::property::prctl::And<Type> const& formula) const {
 		storm::storage::BitVector* result = formula.getLeft().check(*this);
@@ -231,7 +241,7 @@ public:
 	 * Checks the given formula that is a logical "or" of two formulae.
 	 *
 	 * @param formula The formula to check.
-	 * @returns The set of states satisfying the formula represented by a bit vector.
+	 * @return The set of states satisfying the formula represented by a bit vector.
 	 */
 	virtual storm::storage::BitVector* checkOr(storm::property::prctl::Or<Type> const& formula) const {
 		storm::storage::BitVector* result = formula.getLeft().check(*this);
@@ -245,7 +255,7 @@ public:
 	 * Checks the given formula that is a logical "not" of a sub-formula.
 	 *
 	 * @param formula The formula to check.
-	 * @returns The set of states satisfying the formula represented by a bit vector.
+	 * @return The set of states satisfying the formula represented by a bit vector.
 	 */
 	storm::storage::BitVector* checkNot(const storm::property::prctl::Not<Type>& formula) const {
 		storm::storage::BitVector* result = formula.getChild().check(*this);
@@ -258,7 +268,7 @@ public:
 	 * Checks the given formula that is a P operator over a path formula featuring a value bound.
 	 *
 	 * @param formula The formula to check.
-	 * @returns The set of states satisfying the formula represented by a bit vector.
+	 * @return The set of states satisfying the formula represented by a bit vector.
 	 */
 	storm::storage::BitVector* checkProbabilisticBoundOperator(storm::property::prctl::ProbabilisticBoundOperator<Type> const& formula) const {
 		// First, we need to compute the probability for satisfying the path formula for each state.
@@ -284,7 +294,7 @@ public:
 	 * Checks the given formula that is an R operator over a reward formula featuring a value bound.
 	 *
 	 * @param formula The formula to check.
-	 * @returns The set of states satisfying the formula represented by a bit vector.
+	 * @return The set of states satisfying the formula represented by a bit vector.
 	 */
 	storm::storage::BitVector* checkRewardBoundOperator(const storm::property::prctl::RewardBoundOperator<Type>& formula) const {
 		// First, we need to compute the probability for satisfying the path formula for each state.
diff --git a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
index 36f2ff9ee..5ab063ef7 100644
--- a/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseDtmcPrctlModelChecker.h
@@ -88,31 +88,40 @@ public:
 		// First, we need to compute the states that satisfy the sub-formulas of the bounded until-formula.
 		storm::storage::BitVector* leftStates = formula.getLeft().check(*this);
 		storm::storage::BitVector* rightStates = formula.getRight().check(*this);
-
+        std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
+        
         // If we identify the states that have probability 0 of reaching the target states, we can exclude them in the
         // further analysis.
-        storm::storage::BitVector maybeStates = storm::utility::graph::performProbGreater0(this->getModel(), *leftStates, *rightStates, true, formula.getBound());
-        
-        // Now we can eliminate the rows and columns from the original transition probability matrix that have probability 0.
-        storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates);
-        
-        // Compute the new set of target states in the reduced system.
-        storm::storage::BitVector rightStatesInReducedSystem = maybeStates % *rightStates;
+        storm::storage::BitVector statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0(this->getModel(), *leftStates, *rightStates, true, formula.getBound());
+        LOG4CPLUS_INFO(logger, "Found " << statesWithProbabilityGreater0.getNumberOfSetBits() << " 'maybe' states.");
         
-		// Make all rows absorbing that satisfy the second sub-formula.
-		submatrix.makeRowsAbsorbing(rightStatesInReducedSystem);
-
-		// Create the vector with which to multiply.
-        std::vector<Type> subresult(maybeStates.getNumberOfSetBits());
-		storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::constGetOne<Type>());
-
-		// Perform the matrix vector multiplication as often as required by the formula bound.
-		this->performMatrixVectorMultiplication(submatrix, subresult, nullptr, formula.getBound());
-
-        // Create result vector and set its values accordingly.
-		std::vector<Type>* result = new std::vector<Type>(this->getModel().getNumberOfStates());
-        storm::utility::vector::setVectorValues(*result, maybeStates, subresult);
-        storm::utility::vector::setVectorValues<Type>(*result, ~maybeStates, storm::utility::constGetZero<Type>());
+        // Check if we already know the result (i.e. probability 0) for all initial states and
+        // don't compute anything in this case.
+        if (this->getInitialStates().isSubsetOf(~statesWithProbabilityGreater0)) {
+            LOG4CPLUS_INFO(logger, "The probabilities for the initial states were determined in a preprocessing step."
+                           << " No exact probabilities were computed.");
+        } else { // Otherwise, we have have to compute the probabilities.
+            // Now we can eliminate the rows and columns from the original transition probability matrix that have probability 0.
+            storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(statesWithProbabilityGreater0);
+            
+            // Compute the new set of target states in the reduced system.
+            storm::storage::BitVector rightStatesInReducedSystem = statesWithProbabilityGreater0 % *rightStates;
+            
+            // Make all rows absorbing that satisfy the second sub-formula.
+            submatrix.makeRowsAbsorbing(rightStatesInReducedSystem);
+            
+            // Create the vector with which to multiply.
+            std::vector<Type> subresult(statesWithProbabilityGreater0.getNumberOfSetBits());
+            storm::utility::vector::setVectorValues(subresult, rightStatesInReducedSystem, storm::utility::constGetOne<Type>());
+            
+            // Perform the matrix vector multiplication as often as required by the formula bound.
+            this->performMatrixVectorMultiplication(submatrix, subresult, nullptr, formula.getBound());
+            
+            // Create result vector and set its values accordingly.
+
+            storm::utility::vector::setVectorValues(*result, statesWithProbabilityGreater0, subresult);
+            storm::utility::vector::setVectorValues<Type>(*result, ~statesWithProbabilityGreater0, storm::utility::constGetZero<Type>());
+        }
         
 		// Delete obsolete intermediates and return result.
 		delete leftStates;
@@ -143,7 +152,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;
@@ -231,42 +240,49 @@ public:
 		delete rightStates;
 
 		// Perform some logging.
+        storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
 		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 && !qualitative) {
-			// 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();
-
-			// 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(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 = this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, statesWithProbability1);
-
-			// Now solve the created system of linear equations.
-			this->solveEquationSystem(submatrix, x, b);
-
-			// Set values of resulting vector according to result.
-			storm::utility::vector::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.
-			storm::utility::vector::setVectorValues<Type>(*result, maybeStates, Type(0.5));
-		}
+        // Check whether we need to compute exact probabilities for some states.
+        if (this->getInitialStates().isDisjointFrom(maybeStates) || qualitative) {
+            if (qualitative) {
+                LOG4CPLUS_INFO(logger, "The formula was checked qualitatively. No exact probabilities were computed.");
+            } else {
+                LOG4CPLUS_INFO(logger, "The probabilities for the initial states were determined in a preprocessing step."
+                               << " No exact probabilities were computed.");
+            }
+            // Set the values for all maybe-states to 0.5 to indicate that their probability values
+            // are neither 0 nor 1.
+            storm::utility::vector::setVectorValues<Type>(*result, maybeStates, Type(0.5));
+        } else {
+            // Otherwise, we have have to compute the probabilities.
+            
+            // First, 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();
+            
+            // 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(maybeStates.getNumberOfSetBits(), 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 = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(maybeStates, statesWithProbability1);
+            
+            // Now solve the created system of linear equations.
+            this->solveEquationSystem(submatrix, x, b);
+            
+            // Set values of resulting vector according to result.
+            storm::utility::vector::setVectorValues<Type>(*result, maybeStates, x);
+        }
 
 		// Set values of resulting vector that are known exactly.
 		storm::utility::vector::setVectorValues<Type>(*result, statesWithProbability0, storm::utility::constGetZero<Type>());
@@ -294,10 +310,10 @@ public:
 		}
 
 		// Initialize result to state rewards of the model.
-		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
+		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;
@@ -324,24 +340,24 @@ public:
 		// Compute the reward vector to add in each step based on the available reward models.
 		std::vector<Type> totalRewardVector;
 		if (this->getModel().hasTransitionRewards()) {
-			totalRewardVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
+			totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
 			if (this->getModel().hasStateRewards()) {
-				gmm::add(*this->getModel().getStateRewardVector(), totalRewardVector);
+                storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector());
 			}
 		} else {
-			totalRewardVector = std::vector<Type>(*this->getModel().getStateRewardVector());
+			totalRewardVector = std::vector<Type>(this->getModel().getStateRewardVector());
 		}
 
 		// Initialize result to either the state rewards of the model or the null vector.
 		std::vector<Type>* result = nullptr;
 		if (this->getModel().hasStateRewards()) {
-			result = new std::vector<Type>(*this->getModel().getStateRewardVector());
+			result = new std::vector<Type>(this->getModel().getStateRewardVector());
 		} else {
 			result = new std::vector<Type>(this->getModel().getNumberOfStates());
 		}
 
 		// Perform the actual matrix-vector multiplication as long as the bound of the formula is met.
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, &totalRewardVector, formula.getBound());
+		this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, &totalRewardVector, formula.getBound());
 
 		// Return result.
 		return result;
@@ -372,56 +388,67 @@ public:
 		storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true);
 		storm::storage::BitVector infinityStates = storm::utility::graph::performProb1(this->getModel(), trueStates, *targetStates);
 		infinityStates.complement();
+        LOG4CPLUS_INFO(logger, "Found " << infinityStates.getNumberOfSetBits() << " 'infinity' 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.
 		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();
-
-			// 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>());
-
-			// Prepare the right-hand side of the equation system.
-			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::vector::selectVectorValues(b, maybeStates, 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(maybeStatesSetBitCount);
-					storm::utility::vector::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::vector::selectVectorValues(b, maybeStates, *this->getModel().getStateRewardVector());
-			}
-
-			// Now solve the resulting equation system.
-			this->solveEquationSystem(submatrix, x, b);
-
-			// Set values of resulting vector according to result.
-			storm::utility::vector::setVectorValues<Type>(*result, maybeStates, x);
-		}
+        
+        // Check whether we need to compute exact rewards for some states.
+        if (this->getInitialStates().isDisjointFrom(maybeStates)) {
+            LOG4CPLUS_INFO(logger, "The rewards for the initial states were determined in a preprocessing step."
+                            << " No exact rewards were computed.");
+            // Set the values for all maybe-states to 1 to indicate that their reward values
+            // are neither 0 nor infinity.
+            storm::utility::vector::setVectorValues<Type>(*result, maybeStates, storm::utility::constGetOne<Type>());
+        } else {
+            // In this case we have to compute the reward values for the remaining states.
+            const int maybeStatesSetBitCount = maybeStates.getNumberOfSetBits();
+
+            // 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();
+            
+            // 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>());
+            
+            // Prepare the right-hand side of the equation system.
+            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::vector::selectVectorValues(b, maybeStates, 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(maybeStatesSetBitCount);
+                    storm::utility::vector::selectVectorValues(subStateRewards, maybeStates, this->getModel().getStateRewardVector());
+                    storm::utility::vector::addVectorsInPlace(b, 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::vector::selectVectorValues(b, maybeStates, this->getModel().getStateRewardVector());
+            }
+            
+            // Now solve the resulting equation system.
+            this->solveEquationSystem(submatrix, x, b);
+            
+            // Set values of resulting vector according to result.
+            storm::utility::vector::setVectorValues<Type>(*result, maybeStates, x);
+        }
 
 		// Set values of resulting vector that are known exactly.
 		storm::utility::vector::setVectorValues(*result, *targetStates, storm::utility::constGetZero<Type>());
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
index 81c9b45e1..3611a70a2 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
@@ -104,7 +104,7 @@ public:
 		}
         
         // Now we can eliminate the rows and columns from the original transition probability matrix that have probability 0.
-        storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix()->getSubmatrix(maybeStates, *this->getModel().getNondeterministicChoiceIndices());
+        storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices());
 
         // Get the "new" nondeterministic choice indices for the submatrix.
         std::vector<uint_fast64_t> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
@@ -154,7 +154,7 @@ public:
 		// Delete obsolete sub-result.
 		delete nextStates;
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, *this->getModel().getNondeterministicChoiceIndices());
+		this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices());
 
 		// Return result.
 		return result;
@@ -259,7 +259,7 @@ public:
 		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());
+			storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices());
 
 			// Get the "new" nondeterministic choice indices for the submatrix.
 			std::vector<uint_fast64_t> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
@@ -269,7 +269,7 @@ public:
 
 			// 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 = this->getModel().getTransitionMatrix()->getConstrainedRowSumVector(maybeStates, *this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, submatrix.getRowCount());
+			std::vector<Type> b = this->getModel().getTransitionMatrix().getConstrainedRowSumVector(maybeStates, this->getModel().getNondeterministicChoiceIndices(), statesWithProbability1, submatrix.getRowCount());
 
 			// Solve the corresponding system of equations.
 			this->solveEquationSystem(submatrix, x, b, subNondeterministicChoiceIndices);
@@ -304,9 +304,9 @@ public:
 		}
 
 		// Initialize result to state rewards of the model.
-		std::vector<Type>* result = new std::vector<Type>(*this->getModel().getStateRewardVector());
+		std::vector<Type>* result = new std::vector<Type>(this->getModel().getStateRewardVector());
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, *this->getModel().getNondeterministicChoiceIndices(), nullptr, formula.getBound());
+		this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices(), nullptr, formula.getBound());
 
 		// Return result.
 		return result;
@@ -333,23 +333,23 @@ public:
 		// Compute the reward vector to add in each step based on the available reward models.
 		std::vector<Type> totalRewardVector;
 		if (this->getModel().hasTransitionRewards()) {
-			totalRewardVector = this->getModel().getTransitionMatrix()->getPointwiseProductRowSumVector(*this->getModel().getTransitionRewardMatrix());
+			totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
 			if (this->getModel().hasStateRewards()) {
-				gmm::add(*this->getModel().getStateRewardVector(), totalRewardVector);
+                storm::utility::vector::addVectorsInPlace(totalRewardVector, this->getModel().getStateRewardVector());
 			}
 		} else {
-			totalRewardVector = std::vector<Type>(*this->getModel().getStateRewardVector());
+			totalRewardVector = std::vector<Type>(this->getModel().getStateRewardVector());
 		}
 
 		// Initialize result to either the state rewards of the model or the null vector.
 		std::vector<Type>* result = nullptr;
 		if (this->getModel().hasStateRewards()) {
-			result = new std::vector<Type>(*this->getModel().getStateRewardVector());
+			result = new std::vector<Type>(this->getModel().getStateRewardVector());
 		} else {
 			result = new std::vector<Type>(this->getModel().getNumberOfStates());
 		}
 
-		this->performMatrixVectorMultiplication(*this->getModel().getTransitionMatrix(), *result, *this->getModel().getNondeterministicChoiceIndices(), &totalRewardVector, formula.getBound());
+		this->performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), *result, this->getModel().getNondeterministicChoiceIndices(), &totalRewardVector, formula.getBound());
 
 		// Delete temporary variables and return result.
 		return result;
@@ -399,7 +399,7 @@ public:
 		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());
+			storm::storage::SparseMatrix<Type> submatrix = this->getModel().getTransitionMatrix().getSubmatrix(maybeStates, this->getModel().getNondeterministicChoiceIndices());
 
 			// Get the "new" nondeterministic choice indices for the submatrix.
 			std::vector<uint_fast64_t> subNondeterministicChoiceIndices = this->computeNondeterministicChoiceIndicesForConstraint(maybeStates);
@@ -415,8 +415,8 @@ public:
 				// 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::vector::selectVectorValues(b, maybeStates, *this->getModel().getNondeterministicChoiceIndices(), pointwiseProductRowSumVector);
+				std::vector<Type> pointwiseProductRowSumVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
+				storm::utility::vector::selectVectorValues(b, maybeStates, this->getModel().getNondeterministicChoiceIndices(), pointwiseProductRowSumVector);
 
 				if (this->getModel().hasStateRewards()) {
 					// If a state-based reward model is also available, we need to add this vector
@@ -424,15 +424,15 @@ public:
 					// that we still consider (i.e. maybeStates), we need to extract these values
 					// first.
 					std::vector<Type> subStateRewards(b.size());
-					storm::utility::vector::selectVectorValuesRepeatedly(subStateRewards, maybeStates, *this->getModel().getNondeterministicChoiceIndices(), *this->getModel().getStateRewardVector());
-					gmm::add(subStateRewards, b);
+					storm::utility::vector::selectVectorValuesRepeatedly(subStateRewards, maybeStates, this->getModel().getNondeterministicChoiceIndices(), this->getModel().getStateRewardVector());
+					storm::utility::vector::addVectorsInPlace(b, 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::vector::selectVectorValuesRepeatedly(b, maybeStates, *this->getModel().getNondeterministicChoiceIndices(), *this->getModel().getStateRewardVector());
+				storm::utility::vector::selectVectorValuesRepeatedly(b, maybeStates, this->getModel().getNondeterministicChoiceIndices(), this->getModel().getStateRewardVector());
 			}
 
 			// Solve the corresponding system of equations.
@@ -573,7 +573,7 @@ private:
 	 */
 	std::vector<uint_fast64_t> computeNondeterministicChoiceIndicesForConstraint(storm::storage::BitVector const& constraint) const {
 		// First, get a reference to the full nondeterministic choice indices.
-		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = *this->getModel().getNondeterministicChoiceIndices();
+		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = this->getModel().getNondeterministicChoiceIndices();
 
 		// Reserve the known amount of slots for the resulting vector.
 		std::vector<uint_fast64_t> subNondeterministicChoiceIndices(constraint.getNumberOfSetBits() + 1);
diff --git a/src/models/AbstractModel.h b/src/models/AbstractModel.h
index 5d99b0f8f..e632e7879 100644
--- a/src/models/AbstractModel.h
+++ b/src/models/AbstractModel.h
@@ -208,7 +208,7 @@ class AbstractModel: public std::enable_shared_from_this<AbstractModel<T>> {
 		 * @return The size of the state space of the model.
 		 */
 		virtual uint_fast64_t getNumberOfStates() const {
-			return this->getTransitionMatrix()->getColumnCount();
+			return this->getTransitionMatrix().getColumnCount();
 		}
 
 		/*!
@@ -216,7 +216,7 @@ class AbstractModel: public std::enable_shared_from_this<AbstractModel<T>> {
 		 * @return The number of (non-zero) transitions of the model.
 		 */
 		virtual uint_fast64_t getNumberOfTransitions() const {
-			return this->getTransitionMatrix()->getNonZeroEntryCount();
+			return this->getTransitionMatrix().getNonZeroEntryCount();
 		}
 
 		/*!
@@ -245,24 +245,24 @@ class AbstractModel: public std::enable_shared_from_this<AbstractModel<T>> {
 		 * @return A pointer to the matrix representing the transition probability
 		 * function.
 		 */
-		std::shared_ptr<storm::storage::SparseMatrix<T>> getTransitionMatrix() const {
-			return transitionMatrix;
+		storm::storage::SparseMatrix<T> const& getTransitionMatrix() const {
+			return *transitionMatrix;
 		}
 
 		/*!
 		 * Returns a pointer to the matrix representing the transition rewards.
 		 * @return A pointer to the matrix representing the transition rewards.
 		 */
-		std::shared_ptr<storm::storage::SparseMatrix<T>> getTransitionRewardMatrix() const {
-			return transitionRewardMatrix;
+		storm::storage::SparseMatrix<T> const& getTransitionRewardMatrix() const {
+			return *transitionRewardMatrix;
 		}
 
 		/*!
 		 * Returns a pointer to the vector representing the state rewards.
 		 * @return A pointer to the vector representing the state rewards.
 		 */
-		std::shared_ptr<std::vector<T>> getStateRewardVector() const {
-			return stateRewardVector;
+		std::vector<T> const& getStateRewardVector() const {
+			return *stateRewardVector;
 		}
 
 		/*!
@@ -277,8 +277,8 @@ class AbstractModel: public std::enable_shared_from_this<AbstractModel<T>> {
 		 * Returns the state labeling associated with this model.
 		 * @return The state labeling associated with this model.
 		 */
-		std::shared_ptr<storm::models::AtomicPropositionsLabeling> getStateLabeling() const {
-			return stateLabeling;
+		storm::models::AtomicPropositionsLabeling const& getStateLabeling() const {
+			return *stateLabeling;
 		}
 
 		/*!
@@ -323,7 +323,7 @@ class AbstractModel: public std::enable_shared_from_this<AbstractModel<T>> {
 			out << "Model type: \t\t" << this->getType() << std::endl;
 			out << "States: \t\t" << this->getNumberOfStates() << std::endl;
 			out << "Transitions: \t\t" << this->getNumberOfTransitions() << std::endl;
-			this->getStateLabeling()->printAtomicPropositionsInformationToStream(out);
+			this->getStateLabeling().printAtomicPropositionsInformationToStream(out);
 			out << "Size in memory: \t"
 				<< (this->getSizeInMemory())/1024 << " kbytes" << std::endl;
 			out << "-------------------------------------------------------------- "
diff --git a/src/models/AbstractNondeterministicModel.h b/src/models/AbstractNondeterministicModel.h
index b64359d2f..0aacc698b 100644
--- a/src/models/AbstractNondeterministicModel.h
+++ b/src/models/AbstractNondeterministicModel.h
@@ -73,8 +73,8 @@ class AbstractNondeterministicModel: public AbstractModel<T> {
 		 * @param the vector indicating which matrix rows represent non-deterministic choices
 		 * of a certain state.
 		 */
-		std::shared_ptr<std::vector<uint_fast64_t>> getNondeterministicChoiceIndices() const {
-			return nondeterministicChoiceIndices;
+		std::vector<uint_fast64_t> const& getNondeterministicChoiceIndices() const {
+			return *nondeterministicChoiceIndices;
 		}
     
         /*!
diff --git a/src/models/Ctmdp.h b/src/models/Ctmdp.h
index e693c70a1..0d6095ea6 100644
--- a/src/models/Ctmdp.h
+++ b/src/models/Ctmdp.h
@@ -84,8 +84,8 @@ private:
 		// Get the settings object to customize linear solving.
 		storm::settings::Settings* s = storm::settings::instance();
 		double precision = s->get<double>("precision");
-		for (uint_fast64_t row = 0; row < this->getTransitionMatrix()->getRowCount(); row++) {
-			T sum = this->getTransitionMatrix()->getRowSum(row);
+		for (uint_fast64_t row = 0; row < this->getTransitionMatrix().getRowCount(); row++) {
+			T sum = this->getTransitionMatrix().getRowSum(row);
 			if (sum == 0) continue;
 			if (std::abs(sum - 1) > precision) return false;
 		}
diff --git a/src/models/Dtmc.h b/src/models/Dtmc.h
index 41e6bae8a..97399115d 100644
--- a/src/models/Dtmc.h
+++ b/src/models/Dtmc.h
@@ -49,8 +49,8 @@ public:
 			LOG4CPLUS_ERROR(logger, "Probability matrix is invalid.");
 			throw storm::exceptions::InvalidArgumentException() << "Probability matrix is invalid.";
 		}
-		if (this->getTransitionRewardMatrix() != nullptr) {
-			if (!this->getTransitionRewardMatrix()->containsAllPositionsOf(*this->getTransitionMatrix())) {
+		if (this->hasTransitionRewards()) {
+			if (!this->getTransitionRewardMatrix().containsAllPositionsOf(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.";
 			}
@@ -89,14 +89,14 @@ private:
 		storm::settings::Settings* s = storm::settings::instance();
 		double precision = s->get<double>("precision");
 
-		if (this->getTransitionMatrix()->getRowCount() != this->getTransitionMatrix()->getColumnCount()) {
+		if (this->getTransitionMatrix().getRowCount() != this->getTransitionMatrix().getColumnCount()) {
 			// not square
 			LOG4CPLUS_ERROR(logger, "Probability matrix is not square.");
 			return false;
 		}
 
-		for (uint_fast64_t row = 0; row < this->getTransitionMatrix()->getRowCount(); ++row) {
-			T sum = this->getTransitionMatrix()->getRowSum(row);
+		for (uint_fast64_t row = 0; row < this->getTransitionMatrix().getRowCount(); ++row) {
+			T sum = this->getTransitionMatrix().getRowSum(row);
 			if (sum == 0) {
 				LOG4CPLUS_ERROR(logger, "Row " << row << " has sum 0");
 				return false;
diff --git a/src/models/Mdp.h b/src/models/Mdp.h
index 586de6e5e..b2a8d285f 100644
--- a/src/models/Mdp.h
+++ b/src/models/Mdp.h
@@ -86,8 +86,8 @@ private:
 		// Get the settings object to customize linear solving.
 		storm::settings::Settings* s = storm::settings::instance();
 		double precision = s->get<double>("precision");
-		for (uint_fast64_t row = 0; row < this->getTransitionMatrix()->getRowCount(); row++) {
-			T sum = this->getTransitionMatrix()->getRowSum(row);
+		for (uint_fast64_t row = 0; row < this->getTransitionMatrix().getRowCount(); row++) {
+			T sum = this->getTransitionMatrix().getRowSum(row);
 			if (sum == 0) continue;
 			if (std::abs(sum - 1) > precision)  {
 				return false;
diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h
index d412b2eb7..376555e3e 100644
--- a/src/storage/BitVector.h
+++ b/src/storage/BitVector.h
@@ -420,6 +420,7 @@ public:
 	 * Performs a logical "implies" with the given bit vector. In case the sizes of the bit vectors
 	 * do not match, only the matching portion is considered and the overlapping bits
 	 * are set to 0.
+     *
 	 * @param bv A reference to the bit vector to use for the operation.
 	 * @return A bit vector corresponding to the logical "implies" of the two bit vectors.
 	 */
@@ -439,14 +440,15 @@ public:
 	/*!
 	 * Checks whether all bits that are set in the current bit vector are also set in the given bit
 	 * vector.
+     *
 	 * @param bv A reference to the bit vector whose bits are (possibly) a superset of the bits of
 	 * the current bit vector.
-	 * @returns True iff all bits that are set in the current bit vector are also set in the given bit
+	 * @return True iff all bits that are set in the current bit vector are also set in the given bit
 	 * vector.
 	 */
-	bool isContainedIn(BitVector const& bv) const {
+	bool isSubsetOf(BitVector const& bv) const {
 		for (uint_fast64_t i = 0; i < this->bucketCount; ++i) {
-			if ((this->bucketArray[i] & bv.bucketArray[i]) != bv.bucketArray[i]) {
+			if ((this->bucketArray[i] & bv.bucketArray[i]) != this->bucketArray[i]) {
 				return false;
 			}
 		}
@@ -456,6 +458,7 @@ public:
 	/*!
 	 * Checks whether none of the bits that are set in the current bit vector are also set in the
 	 * given bit vector.
+     *
 	 * @param bv A reference to the bit vector whose bits are (possibly) disjoint from the bits in
 	 * the current bit vector.
 	 * @returns True iff none of the bits that are set in the current bit vector are also set in the
@@ -504,9 +507,25 @@ public:
         
 		return result;
 	}
+    
+    /*!
+     * Returns a list containing all indices such that the bits at these indices are set to true
+     * in the bit vector.
+     *
+     * @return A vector of indices of set bits in the bit vector.
+     */
+    std::vector<uint_fast64_t> getSetIndicesList() const {
+        std::vector<uint_fast64_t> result;
+        result.reserve(this->getNumberOfSetBits());
+        for (auto index : *this) {
+			result.push_back(index);
+		}
+        return result;
+    }
 
 	/*!
-	 * Adds all indices of bits set to one to the provided list.
+	 * Adds all indices of bits set to one to the given list.
+     *
 	 * @param list The list to which to append the indices.
 	 */
 	void addSetIndicesToVector(std::vector<uint_fast64_t>& vector) const {
diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h
index 4ee005881..96d56424e 100644
--- a/src/storage/SparseMatrix.h
+++ b/src/storage/SparseMatrix.h
@@ -920,7 +920,7 @@ public:
 	 * @returns A vector containing the sum of the elements in each row of the matrix resulting from
 	 * pointwise multiplication of the current matrix with the given matrix.
 	 */
-	std::vector<T> getPointwiseProductRowSumVector(storm::storage::SparseMatrix<T> const& otherMatrix) {
+	std::vector<T> getPointwiseProductRowSumVector(storm::storage::SparseMatrix<T> const& otherMatrix) const {
 		// Prepare result.
 		std::vector<T> result(rowCount, storm::utility::constGetZero<T>());
 
diff --git a/src/utility/graph.h b/src/utility/graph.h
index 8d9f0f024..4d67f35ac 100644
--- a/src/utility/graph.h
+++ b/src/utility/graph.h
@@ -46,10 +46,8 @@ namespace graph {
 		// Add all psi states as the already satisfy the condition.
 		statesWithProbabilityGreater0 |= psiStates;
 
-		// Initialize the stack used for the BFS with the states.
-		std::vector<uint_fast64_t> stack;
-		stack.reserve(model.getNumberOfStates());
-		psiStates.addSetIndicesToVector(stack);
+		// Initialize the stack used for the DFS with the states.
+		std::vector<uint_fast64_t> stack = psiStates.getSetIndicesList();
         
         // Initialize the stack for the step bound, if the number of steps is bounded.
         std::vector<uint_fast64_t> stepStack;
@@ -63,7 +61,7 @@ namespace graph {
             }
         }
 
-		// Perform the actual BFS.
+		// Perform the actual DFS.
         uint_fast64_t currentState, currentStepBound;
 		while (!stack.empty()) {
             currentState = stack.back();
@@ -74,17 +72,17 @@ namespace graph {
                 stepStack.pop_back();
             }
 
-			for (auto it = backwardTransitions.constColumnIteratorBegin(currentState); it != backwardTransitions.constColumnIteratorEnd(currentState); ++it) {
-				if (phiStates.get(*it) && (!statesWithProbabilityGreater0.get(*it) || (useStepBound && remainingSteps[*it] < currentStepBound - 1))) {
-                    // If we don't have a number of maximal steps to take, just add the state to the stack.
+			for (auto predecessorIt = backwardTransitions.constColumnIteratorBegin(currentState), predecessorIte = backwardTransitions.constColumnIteratorEnd(currentState); predecessorIt != predecessorIte; ++predecessorIt) {
+				if (phiStates.get(*predecessorIt) && (!statesWithProbabilityGreater0.get(*predecessorIt) || (useStepBound && remainingSteps[*predecessorIt] < currentStepBound - 1))) {
+                    // If we don't have a bound on the number of steps to take, just add the state to the stack.
                     if (!useStepBound) {
-                        statesWithProbabilityGreater0.set(*it, true);
-                        stack.push_back(*it);
+                        statesWithProbabilityGreater0.set(*predecessorIt, true);
+                        stack.push_back(*predecessorIt);
                     } else if (currentStepBound > 0) {
                         // If there is at least one more step to go, we need to push the state and the new number of steps.
-                        remainingSteps[*it] = currentStepBound - 1;
-                        statesWithProbabilityGreater0.set(*it, true);
-                        stack.push_back(*it);
+                        remainingSteps[*predecessorIt] = currentStepBound - 1;
+                        statesWithProbabilityGreater0.set(*predecessorIt, true);
+                        stack.push_back(*predecessorIt);
                         stepStack.push_back(currentStepBound - 1);
                     }
 				}
@@ -174,17 +172,11 @@ namespace graph {
         // Prepare resulting bit vector.
         storm::storage::BitVector statesWithProbabilityGreater0(model.getNumberOfStates());
         
-		// Get some temporaries for convenience.
-		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix = model.getTransitionMatrix();
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
-        
 		// Add all psi states as the already satisfy the condition.
 		statesWithProbabilityGreater0 |= psiStates;
         
-		// Initialize the stack used for the BFS with the states
-		std::vector<uint_fast64_t> stack;
-		stack.reserve(model.getNumberOfStates());
-		psiStates.addSetIndicesToVector(stack);
+		// Initialize the stack used for the DFS with the states
+		std::vector<uint_fast64_t> stack = psiStates.getSetIndicesList();
         
         // Initialize the stack for the step bound, if the number of steps is bounded.
         std::vector<uint_fast64_t> stepStack;
@@ -198,7 +190,7 @@ namespace graph {
             }
         }
         
-		// Perform the actual BFS.
+		// Perform the actual DFS.
         uint_fast64_t currentState, currentStepBound;
 		while (!stack.empty()) {
 			currentState = stack.back();
@@ -209,17 +201,17 @@ namespace graph {
                 stepStack.pop_back();
             }
             
-			for (auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
-                if (phiStates.get(*it) && (!statesWithProbabilityGreater0.get(*it) || (useStepBound && remainingSteps[*it] < currentStepBound - 1))) {
+			for (auto predecessorIt = backwardTransitions.constColumnIteratorBegin(currentState), predecessorIte = backwardTransitions.constColumnIteratorEnd(currentState); predecessorIt != predecessorIte; ++predecessorIt) {
+                if (phiStates.get(*predecessorIt) && (!statesWithProbabilityGreater0.get(*predecessorIt) || (useStepBound && remainingSteps[*predecessorIt] < currentStepBound - 1))) {
                     // If we don't have a bound on the number of steps to take, just add the state to the stack.
                     if (!useStepBound) {
-                        statesWithProbabilityGreater0.set(*it, true);
-                        stack.push_back(*it);
+                        statesWithProbabilityGreater0.set(*predecessorIt, true);
+                        stack.push_back(*predecessorIt);
                     } else if (currentStepBound > 0) {
                         // If there is at least one more step to go, we need to push the state and the new number of steps.
-                        remainingSteps[*it] = currentStepBound - 1;
-                        statesWithProbabilityGreater0.set(*it, true);
-                        stack.push_back(*it);
+                        remainingSteps[*predecessorIt] = currentStepBound - 1;
+                        statesWithProbabilityGreater0.set(*predecessorIt, true);
+                        stack.push_back(*predecessorIt);
                         stepStack.push_back(currentStepBound - 1);
                     }
                 }
@@ -265,11 +257,11 @@ namespace graph {
 	template <typename T>
 	storm::storage::BitVector performProb1E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<bool> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
         // Get some temporaries for convenience.
-		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix = model.getTransitionMatrix();
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
+        storm::storage::SparseMatrix<T> const& transitionMatrix = model.getTransitionMatrix();
+		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
 
+        // Initialize the environment for the iterative algorithm.
 		storm::storage::BitVector currentStates(model.getNumberOfStates(), true);
-
 		std::vector<uint_fast64_t> stack;
 		stack.reserve(model.getNumberOfStates());
 
@@ -285,14 +277,14 @@ namespace graph {
 				currentState = stack.back();
 				stack.pop_back();
 
-				for(auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
-					if (phiStates.get(*it) && !nextStates.get(*it)) {
+				for(auto predecessorIt = backwardTransitions.constColumnIteratorBegin(currentState), predecessorIte = backwardTransitions.constColumnIteratorEnd(currentState); predecessorIt != predecessorIte; ++predecessorIt) {
+					if (phiStates.get(*predecessorIt) && !nextStates.get(*predecessorIt)) {
 						// Check whether the predecessor has only successors in the current state set for one of the
 						// nondeterminstic choices.
-						for (uint_fast64_t row = (*nondeterministicChoiceIndices)[*it]; row < (*nondeterministicChoiceIndices)[*it + 1]; ++row) {
+						for (auto row = nondeterministicChoiceIndices[*predecessorIt]; row < nondeterministicChoiceIndices[*predecessorIt + 1]; ++row) {
 							bool allSuccessorsInCurrentStates = true;
-							for (auto colIt = transitionMatrix->constColumnIteratorBegin(row); colIt != transitionMatrix->constColumnIteratorEnd(row); ++colIt) {
-								if (!currentStates.get(*colIt)) {
+							for (auto targetIt = transitionMatrix.constColumnIteratorBegin(row), targetIte = transitionMatrix.constColumnIteratorEnd(row); targetIt != targetIte; ++targetIt) {
+								if (!currentStates.get(*targetIt)) {
 									allSuccessorsInCurrentStates = false;
 									break;
 								}
@@ -302,8 +294,8 @@ namespace graph {
 							// add it to the set of states for the next iteration and perform a backward search from
 							// that state.
 							if (allSuccessorsInCurrentStates) {
-								nextStates.set(*it, true);
-								stack.push_back(*it);
+								nextStates.set(*predecessorIt, true);
+								stack.push_back(*predecessorIt);
 								break;
 							}
 						}
@@ -364,16 +356,14 @@ namespace graph {
         storm::storage::BitVector statesWithProbabilityGreater0(model.getNumberOfStates());
         
 		// Get some temporaries for convenience.
-		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix = model.getTransitionMatrix();
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
+        storm::storage::SparseMatrix<T> const& transitionMatrix = model.getTransitionMatrix();
+		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
         
 		// Add all psi states as the already satisfy the condition.
 		statesWithProbabilityGreater0 |= psiStates;
         
-		// Initialize the stack used for the BFS with the states
-		std::vector<uint_fast64_t> stack;
-		stack.reserve(model.getNumberOfStates());
-		psiStates.addSetIndicesToVector(stack);
+		// Initialize the stack used for the DFS with the states
+		std::vector<uint_fast64_t> stack = psiStates.getSetIndicesList();
         
         // Initialize the stack for the step bound, if the number of steps is bounded.
         std::vector<uint_fast64_t> stepStack;
@@ -387,7 +377,7 @@ namespace graph {
             }
         }
         
-		// Perform the actual BFS.
+		// Perform the actual DFS.
         uint_fast64_t currentState, currentStepBound;
 		while(!stack.empty()) {
 			currentState = stack.back();
@@ -398,15 +388,15 @@ namespace graph {
                 stepStack.pop_back();
             }
             
-			for(auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
-                if (phiStates.get(*it) && (!statesWithProbabilityGreater0.get(*it) || (useStepBound && remainingSteps[*it] < currentStepBound - 1))) {
+			for(auto predecessorIt = backwardTransitions.constColumnIteratorBegin(currentState), predecessorIte = backwardTransitions.constColumnIteratorEnd(currentState); predecessorIt != predecessorIte; ++predecessorIt) {
+                if (phiStates.get(*predecessorIt) && (!statesWithProbabilityGreater0.get(*predecessorIt) || (useStepBound && remainingSteps[*predecessorIt] < currentStepBound - 1))) {
                     // Check whether the predecessor has at least one successor in the current state set for every
                     // nondeterministic choice.
                     bool addToStatesWithProbabilityGreater0 = true;
-                    for (auto rowIt = nondeterministicChoiceIndices->begin() + *it; rowIt != nondeterministicChoiceIndices->begin() + *it + 1; ++rowIt) {
+                    for (auto row = nondeterministicChoiceIndices[*predecessorIt]; row < nondeterministicChoiceIndices[*predecessorIt + 1]; ++row) {
                         bool hasAtLeastOneSuccessorWithProbabilityGreater0 = false;
-                        for (auto colIt = transitionMatrix->constColumnIteratorBegin(*rowIt); colIt != transitionMatrix->constColumnIteratorEnd(*rowIt); ++colIt) {
-                            if (statesWithProbabilityGreater0.get(*colIt)) {
+                        for (auto successorIt = transitionMatrix.constColumnIteratorBegin(row), successorIte = transitionMatrix.constColumnIteratorEnd(row); successorIt != successorIte; ++successorIt) {
+                            if (statesWithProbabilityGreater0.get(*successorIt)) {
                                 hasAtLeastOneSuccessorWithProbabilityGreater0 = true;
                                 break;
                             }
@@ -422,13 +412,13 @@ namespace graph {
                     if (addToStatesWithProbabilityGreater0) {
                         // If we don't have a bound on the number of steps to take, just add the state to the stack.
                         if (!useStepBound) {
-                            statesWithProbabilityGreater0.set(*it, true);
-                            stack.push_back(*it);
+                            statesWithProbabilityGreater0.set(*predecessorIt, true);
+                            stack.push_back(*predecessorIt);
                         } else if (currentStepBound > 0) {
                             // If there is at least one more step to go, we need to push the state and the new number of steps.
-                            remainingSteps[*it] = currentStepBound - 1;
-                            statesWithProbabilityGreater0.set(*it, true);
-                            stack.push_back(*it);
+                            remainingSteps[*predecessorIt] = currentStepBound - 1;
+                            statesWithProbabilityGreater0.set(*predecessorIt, true);
+                            stack.push_back(*predecessorIt);
                             stepStack.push_back(currentStepBound - 1);
                         }
                     }
@@ -473,8 +463,8 @@ namespace graph {
 	template <typename T>
 	storm::storage::BitVector performProb1A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<bool> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
 		// Get some temporaries for convenience.
-		std::shared_ptr<storm::storage::SparseMatrix<T>> transitionMatrix = model.getTransitionMatrix();
-		std::shared_ptr<std::vector<uint_fast64_t>> nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
+        storm::storage::SparseMatrix<T> const& transitionMatrix = model.getTransitionMatrix();
+		std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
 
         // Initialize the environment for the iterative algorithm.
 		storm::storage::BitVector currentStates(model.getNumberOfStates(), true);
@@ -493,14 +483,14 @@ namespace graph {
 				currentState = stack.back();
 				stack.pop_back();
 
-				for(auto it = backwardTransitions.constColumnIteratorBegin(currentState), ite = backwardTransitions.constColumnIteratorEnd(currentState); it != ite; ++it) {
-					if (phiStates.get(*it) && !nextStates.get(*it)) {
+				for(auto predecessorIt = backwardTransitions.constColumnIteratorBegin(currentState), predecessorIte = backwardTransitions.constColumnIteratorEnd(currentState); predecessorIt != predecessorIte; ++predecessorIt) {
+					if (phiStates.get(*predecessorIt) && !nextStates.get(*predecessorIt)) {
 						// Check whether the predecessor has only successors in the current state set for all of the
 						// nondeterminstic choices.
 						bool allSuccessorsInCurrentStatesForAllChoices = true;
-						for (uint_fast64_t row = (*nondeterministicChoiceIndices)[*it]; row < (*nondeterministicChoiceIndices)[*it + 1]; ++row) {
-							for (auto colIt = transitionMatrix->constColumnIteratorBegin(row); colIt != transitionMatrix->constColumnIteratorEnd(row); ++colIt) {
-								if (!currentStates.get(*colIt)) {
+						for (auto row = nondeterministicChoiceIndices[*predecessorIt]; row < nondeterministicChoiceIndices[*predecessorIt + 1]; ++row) {
+							for (auto successorIt = transitionMatrix.constColumnIteratorBegin(row), successorIte = transitionMatrix.constColumnIteratorEnd(row); successorIt != successorIte; ++successorIt) {
+								if (!currentStates.get(*successorIt)) {
 									allSuccessorsInCurrentStatesForAllChoices = false;
 									goto afterCheckLoop;
 								}
@@ -512,8 +502,8 @@ namespace graph {
 						// add it to the set of states for the next iteration and perform a backward search from
 						// that state.
 						if (allSuccessorsInCurrentStatesForAllChoices) {
-							nextStates.set(*it, true);
-							stack.push_back(*it);
+							nextStates.set(*predecessorIt, true);
+							stack.push_back(*predecessorIt);
 						}
 					}
 				}
@@ -584,7 +574,7 @@ namespace graph {
     recursionStepForward:
 		while (!recursionStateStack.empty()) {
 			uint_fast64_t currentState = recursionStateStack.back();
-			typename storm::storage::SparseMatrix<T>::ConstIndexIterator currentIt = recursionIteratorStack.back();
+			typename storm::storage::SparseMatrix<T>::ConstIndexIterator successorIt = recursionIteratorStack.back();
             
 			// Perform the treatment of newly discovered state as defined by Tarjan's algorithm
 			visitedStates.set(currentState, true);
@@ -595,29 +585,29 @@ namespace graph {
 			tarjanStackStates.set(currentState, true);
             
 			// Now, traverse all successors of the current state.
-			for(; currentIt != matrix.constColumnIteratorEnd(currentState); ++currentIt) {
+			for(; successorIt != matrix.constColumnIteratorEnd(currentState); ++successorIt) {
 				// If we have not visited the successor already, we need to perform the procedure
 				// recursively on the newly found state.
-				if (!visitedStates.get(*currentIt)) {
+				if (!visitedStates.get(*successorIt)) {
 					// Save current iterator position so we can continue where we left off later.
 					recursionIteratorStack.pop_back();
-					recursionIteratorStack.push_back(currentIt);
+					recursionIteratorStack.push_back(successorIt);
                     
 					// Put unvisited successor on top of our recursion stack and remember that.
-					recursionStateStack.push_back(*currentIt);
-					statesInStack[*currentIt] = true;
+					recursionStateStack.push_back(*successorIt);
+					statesInStack[*successorIt] = true;
                     
 					// Also, put initial value for iterator on corresponding recursion stack.
-					recursionIteratorStack.push_back(matrix.constColumnIteratorBegin(*currentIt));
+					recursionIteratorStack.push_back(matrix.constColumnIteratorBegin(*successorIt));
                     
 					// Perform the actual recursion step in an iterative way.
 					goto recursionStepForward;
                     
                 recursionStepBackward:
-					lowlinks[currentState] = std::min(lowlinks[currentState], lowlinks[*currentIt]);
-				} else if (tarjanStackStates.get(*currentIt)) {
+					lowlinks[currentState] = std::min(lowlinks[currentState], lowlinks[*successorIt]);
+				} else if (tarjanStackStates.get(*successorIt)) {
 					// Update the lowlink of the current state.
-					lowlinks[currentState] = std::min(lowlinks[currentState], stateIndices[*currentIt]);
+					lowlinks[currentState] = std::min(lowlinks[currentState], stateIndices[*successorIt]);
 				}
 			}
             
@@ -648,7 +638,7 @@ namespace graph {
 			// original recursive call.
 			if (recursionStateStack.size() > 0) {
 				currentState = recursionStateStack.back();
-				currentIt = recursionIteratorStack.back();
+				successorIt = recursionIteratorStack.back();
                 
 				goto recursionStepBackward;
 			}
@@ -680,7 +670,7 @@ namespace graph {
         uint_fast64_t currentIndex = 0;
         for (uint_fast64_t state = 0; state < numberOfStates; ++state) {
             if (!visitedStates.get(state)) {
-                performSccDecompositionHelper(state, currentIndex, stateIndices, lowlinks, tarjanStack, tarjanStackStates, visitedStates, *model.getTransitionMatrix(), scc);
+                performSccDecompositionHelper(state, currentIndex, stateIndices, lowlinks, tarjanStack, tarjanStackStates, visitedStates, model.getTransitionMatrix(), scc);
             }
         }
 
@@ -723,22 +713,22 @@ namespace graph {
 				recursionStepForward:
 				while (!recursionStack.empty()) {
 					uint_fast64_t currentState = recursionStack.back();
-					typename storm::storage::SparseMatrix<T>::ConstIndexIterator currentIt = iteratorRecursionStack.back();
+					typename storm::storage::SparseMatrix<T>::ConstIndexIterator successorIt = iteratorRecursionStack.back();
 
 					visitedStates.set(currentState, true);
 
 					recursionStepBackward:
-					for (; currentIt != matrix.constColumnIteratorEnd(currentState); ++currentIt) {
-						if (!visitedStates.get(*currentIt)) {
+					for (; successorIt != matrix.constColumnIteratorEnd(currentState); ++successorIt) {
+						if (!visitedStates.get(*successorIt)) {
 							// Put unvisited successor on top of our recursion stack and remember that.
-							recursionStack.push_back(*currentIt);
+							recursionStack.push_back(*successorIt);
 
 							// Save current iterator position so we can continue where we left off later.
 							iteratorRecursionStack.pop_back();
-							iteratorRecursionStack.push_back(currentIt + 1);
+							iteratorRecursionStack.push_back(successorIt + 1);
 
 							// Also, put initial value for iterator on corresponding recursion stack.
-							iteratorRecursionStack.push_back(matrix.constColumnIteratorBegin(*currentIt));
+							iteratorRecursionStack.push_back(matrix.constColumnIteratorBegin(*successorIt));
 
 							goto recursionStepForward;
 						}
@@ -756,7 +746,7 @@ namespace graph {
 					// original recursive call.
 					if (recursionStack.size() > 0) {
 						currentState = recursionStack.back();
-						currentIt = iteratorRecursionStack.back();
+						successorIt = iteratorRecursionStack.back();
 
 						goto recursionStepBackward;
 					}
diff --git a/test/functional/parser/ParseMdpTest.cpp b/test/functional/parser/ParseMdpTest.cpp
index dfaa0ac15..a379e82c2 100644
--- a/test/functional/parser/ParseMdpTest.cpp
+++ b/test/functional/parser/ParseMdpTest.cpp
@@ -17,12 +17,12 @@ TEST(ParseMdpTest, parseAndOutput) {
 			STORM_CPP_TESTS_BASE_PATH "/functional/parser/lab_files/pctl_general_input_01.lab"));
 
 	std::shared_ptr<storm::models::Mdp<double>> mdp = mdpParser->getMdp();
-	std::shared_ptr<storm::storage::SparseMatrix<double>> matrix = mdp->getTransitionMatrix();
+	storm::storage::SparseMatrix<double> const& matrix = mdp->getTransitionMatrix();
 
 	ASSERT_EQ(mdp->getNumberOfStates(), (uint_fast64_t)3);
 	ASSERT_EQ(mdp->getNumberOfTransitions(), (uint_fast64_t)11);
-	ASSERT_EQ(matrix->getRowCount(), (uint_fast64_t)(2 * 3));
-	ASSERT_EQ(matrix->getColumnCount(), (uint_fast64_t)3);
+	ASSERT_EQ(matrix.getRowCount(), (uint_fast64_t)(2 * 3));
+	ASSERT_EQ(matrix.getColumnCount(), (uint_fast64_t)3);
 	
 
 	delete mdpParser;