diff --git a/src/counterexamples/MILPMinimalLabelSetGenerator.h b/src/counterexamples/MILPMinimalLabelSetGenerator.h
index 5023c342e..9b3c2b3e2 100644
--- a/src/counterexamples/MILPMinimalLabelSetGenerator.h
+++ b/src/counterexamples/MILPMinimalLabelSetGenerator.h
@@ -90,10 +90,9 @@ namespace storm {
              */
             static struct StateInformation determineRelevantAndProblematicStates(storm::models::Mdp<T> const& labeledMdp, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
                 StateInformation result;
-                storm::storage::SparseMatrix<T> backwardTransitions = labeledMdp.getBackwardTransitions();
-                result.relevantStates = storm::utility::graph::performProbGreater0E(labeledMdp, backwardTransitions, phiStates, psiStates);
+                result.relevantStates = storm::utility::graph::performProbGreater0E(labeledMdp.getTransitionMatrix(), labeledMdp.getNondeterministicChoiceIndices(), labeledMdp.getBackwardTransitions(), phiStates, psiStates);
                 result.relevantStates &= ~psiStates;
-                result.problematicStates = storm::utility::graph::performProbGreater0E(labeledMdp, backwardTransitions, phiStates, psiStates);
+                result.problematicStates = storm::utility::graph::performProbGreater0E(labeledMdp.getTransitionMatrix(), labeledMdp.getNondeterministicChoiceIndices(), labeledMdp.getBackwardTransitions(), phiStates, psiStates);
                 result.problematicStates.complement();
                 result.problematicStates &= result.relevantStates;
                 LOG4CPLUS_DEBUG(logger, "Found " << phiStates.getNumberOfSetBits() << " filter states.");
diff --git a/src/counterexamples/PathBasedSubsystemGenerator.h b/src/counterexamples/PathBasedSubsystemGenerator.h
index 18355e1c4..976da3120 100644
--- a/src/counterexamples/PathBasedSubsystemGenerator.h
+++ b/src/counterexamples/PathBasedSubsystemGenerator.h
@@ -570,7 +570,7 @@ public:
 		storm::property::prctl::Until<T> const* until = dynamic_cast<storm::property::prctl::Until<T> const*>(pathFormulaPtr);
 		if(eventually != nullptr) {
 			targetStates = eventually->getChild().check(modelCheck);
-			allowedStates = storm::storage::BitVector(targetStates.getSize(), true);
+			allowedStates = storm::storage::BitVector(targetStates.size(), true);
 		}
 		else if(globally != nullptr){
 			//eventually reaching a state without property visiting only states with property
diff --git a/src/counterexamples/SMTMinimalCommandSetGenerator.h b/src/counterexamples/SMTMinimalCommandSetGenerator.h
index a4041a070..bd8d45499 100644
--- a/src/counterexamples/SMTMinimalCommandSetGenerator.h
+++ b/src/counterexamples/SMTMinimalCommandSetGenerator.h
@@ -103,7 +103,7 @@ namespace storm {
                 // Compute all relevant states, i.e. states for which there exists a scheduler that has a non-zero
                 // probabilitiy of satisfying phi until psi.
                 storm::storage::SparseMatrix<T> backwardTransitions = labeledMdp.getBackwardTransitions();
-                relevancyInformation.relevantStates = storm::utility::graph::performProbGreater0E(labeledMdp, backwardTransitions, phiStates, psiStates);
+                relevancyInformation.relevantStates = storm::utility::graph::performProbGreater0E(labeledMdp.getTransitionMatrix(), labeledMdp.getNondeterministicChoiceIndices(), backwardTransitions, phiStates, psiStates);
                 relevancyInformation.relevantStates &= ~psiStates;
 
                 LOG4CPLUS_DEBUG(logger, "Found " << relevancyInformation.relevantStates.getNumberOfSetBits() << " relevant states.");
@@ -1427,7 +1427,7 @@ namespace storm {
                 }
                 
                 storm::storage::BitVector unreachableRelevantStates = ~reachableStates & relevancyInformation.relevantStates;
-                storm::storage::BitVector statesThatCanReachTargetStates = storm::utility::graph::performProbGreater0E(subMdp, subMdp.getBackwardTransitions(), phiStates, psiStates);
+                storm::storage::BitVector statesThatCanReachTargetStates = storm::utility::graph::performProbGreater0E(subMdp.getTransitionMatrix(), subMdp.getNondeterministicChoiceIndices(), subMdp.getBackwardTransitions(), phiStates, psiStates);
                 std::vector<storm::storage::VectorSet<uint_fast64_t>> guaranteedLabelSets = storm::utility::counterexamples::getGuaranteedLabelSets(originalMdp, statesThatCanReachTargetStates, relevancyInformation.relevantLabels);
                 
                 LOG4CPLUS_DEBUG(logger, "Found " << reachableLabels.size() << " reachable labels and " << reachableStates.getNumberOfSetBits() << " reachable states.");
@@ -1547,7 +1547,7 @@ namespace storm {
                 LOG4CPLUS_DEBUG(logger, "Successfully determined reachable state space.");
                 
                 storm::storage::BitVector unreachableRelevantStates = ~reachableStates & relevancyInformation.relevantStates;
-                storm::storage::BitVector statesThatCanReachTargetStates = storm::utility::graph::performProbGreater0E(subMdp, subMdp.getBackwardTransitions(), phiStates, psiStates);
+                storm::storage::BitVector statesThatCanReachTargetStates = storm::utility::graph::performProbGreater0E(subMdp.getTransitionMatrix(), subMdp.getNondeterministicChoiceIndices(), subMdp.getBackwardTransitions(), phiStates, psiStates);
                 std::vector<storm::storage::VectorSet<uint_fast64_t>> guaranteedLabelSets = storm::utility::counterexamples::getGuaranteedLabelSets(originalMdp, statesThatCanReachTargetStates, relevancyInformation.relevantLabels);
                 
                 // Search for states for which we could enable another option and possibly improve the reachability probability.
diff --git a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
index 80943c8e0..f050a2286 100644
--- a/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
+++ b/src/modelchecker/csl/SparseMarkovAutomatonCslModelChecker.h
@@ -4,6 +4,7 @@
 #include <stack>
 
 #include "src/modelchecker/csl/AbstractModelChecker.h"
+#include "src/modelchecker/prctl/SparseMdpPrctlModelChecker.h"
 #include "src/models/MarkovAutomaton.h"
 #include "src/storage/BitVector.h"
 #include "src/storage/MaximalEndComponentDecomposition.h"
@@ -20,7 +21,6 @@ namespace storm {
             template<typename ValueType>
             class SparseMarkovAutomatonCslModelChecker : public AbstractModelChecker<ValueType> {
             public:
-                
                 explicit SparseMarkovAutomatonCslModelChecker(storm::models::MarkovAutomaton<ValueType> const& model, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<ValueType>> nondeterministicLinearEquationSolver = storm::utility::solver::getNondeterministicLinearEquationSolver<ValueType>()) : AbstractModelChecker<ValueType>(model), minimumOperatorStack(), nondeterministicLinearEquationSolver(nondeterministicLinearEquationSolver) {
                     // Intentionally left empty.
                 }
@@ -34,7 +34,13 @@ namespace storm {
                 }
                 
                 std::vector<ValueType> checkUntil(storm::property::csl::Until<ValueType> const& formula, bool qualitative) const {
-                    throw storm::exceptions::NotImplementedException() << "Model checking Until formulas on Markov automata is not yet implemented.";
+                    storm::storage::BitVector leftStates = formula.getLeft().check(*this);
+                    storm::storage::BitVector rightStates = formula.getRight().check(*this);
+                    return computeUnboundedUntilProbabilities(minimumOperatorStack.top(), leftStates, rightStates, qualitative).first;
+                }
+                
+                std::pair<std::vector<ValueType>, storm::storage::TotalScheduler> computeUnboundedUntilProbabilities(bool min, storm::storage::BitVector const& leftStates, storm::storage::BitVector const& rightStates, bool qualitative) const {
+                    return storm::modelchecker::prctl::SparseMdpPrctlModelChecker<ValueType>::computeUnboundedUntilProbabilities(min, this->getModel().getTransitionMatrix(), this->getModel().getNondeterministicChoiceIndices(), this->getModel().getBackwardTransitions(), this->getModel().getInitialStates(), leftStates, rightStates, nondeterministicLinearEquationSolver, qualitative);
                 }
                 
                 std::vector<ValueType> checkTimeBoundedUntil(storm::property::csl::TimeBoundedUntil<ValueType> const& formula, bool qualitative) const {
@@ -42,7 +48,8 @@ namespace storm {
                 }
                 
                 std::vector<ValueType> checkTimeBoundedEventually(storm::property::csl::TimeBoundedEventually<ValueType> const& formula, bool qualitative) const {
-                    throw storm::exceptions::NotImplementedException() << "Model checking time-bounded Until formulas on Markov automata is not yet implemented.";
+                    storm::storage::BitVector goalStates = formula.getChild().check(*this);
+                    return this->checkTimeBoundedEventually(this->minimumOperatorStack.top(), goalStates, formula.getLowerBound(), formula.getUpperBound());
                 }
                 
                 std::vector<ValueType> checkGlobally(storm::property::csl::Globally<ValueType> const& formula, bool qualitative) const {
@@ -50,7 +57,8 @@ namespace storm {
                 }
                 
                 std::vector<ValueType> checkEventually(storm::property::csl::Eventually<ValueType> const& formula, bool qualitative) const {
-                    throw storm::exceptions::NotImplementedException() << "Model checking Eventually formulas on Markov automata is not yet implemented.";
+                    storm::storage::BitVector subFormulaStates = formula.getChild().check(*this);
+                    return computeUnboundedUntilProbabilities(minimumOperatorStack.top(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), subFormulaStates, qualitative).first;
                 }
                 
                 std::vector<ValueType> checkNext(storm::property::csl::Next<ValueType> const& formula, bool qualitative) const {
@@ -391,18 +399,13 @@ namespace storm {
                     // Prepare result vector.
                     std::vector<ValueType> result(this->getModel().getNumberOfStates());
                     
-                    std::cout << "res " << result << std::endl;
-                    
                     // Set the values for states not contained in MECs.
-                    std::cout << x << std::endl;
                     storm::utility::vector::setVectorValues(result, statesNotContainedInAnyMec, x);
-                    std::cout << "res " << result << std::endl;
                     
                     // Set the values for all states in MECs.
                     for (auto state : statesInMecs) {
                         result[state] = lraValuesForEndComponents[stateToMecIndexMap[state]];
                     }
-                    std::cout << "res " << result << std::endl;
 
                     return result;
                 }
diff --git a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
index b51bb5565..008ea3606 100644
--- a/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
+++ b/src/modelchecker/prctl/SparseMdpPrctlModelChecker.h
@@ -96,9 +96,9 @@ namespace storm {
                     // Determine the states that have 0 probability of reaching the target states.
                     storm::storage::BitVector statesWithProbabilityGreater0;
                     if (this->minimumOperatorStack.top()) {
-                        statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0A(this->getModel(), this->getModel().getBackwardTransitions(), phiStates, psiStates, true, stepBound);
+                        statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0A(this->getModel().getTransitionMatrix(), this->getModel().getNondeterministicChoiceIndices(), this->getModel().getBackwardTransitions(), phiStates, psiStates, true, stepBound);
                     } else {
-                        statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0E(this->getModel(), this->getModel().getBackwardTransitions(), phiStates, psiStates, true, stepBound);
+                        statesWithProbabilityGreater0 = storm::utility::graph::performProbGreater0E(this->getModel().getTransitionMatrix(), this->getModel().getNondeterministicChoiceIndices(), this->getModel().getBackwardTransitions(), phiStates, psiStates, true, stepBound);
                     }
                     
                     // Check if we already know the result (i.e. probability 0) for all initial states and
@@ -278,29 +278,30 @@ namespace storm {
                  * @return The probabilities for the satisfying phi until psi for each state of the model. If the
                  * qualitative flag is set, exact probabilities might not be computed.
                  */
-                std::pair<std::vector<Type>, storm::storage::TotalScheduler> checkUntil(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const {
+                static std::pair<std::vector<Type>, storm::storage::TotalScheduler> computeUnboundedUntilProbabilities(bool minimize, storm::storage::SparseMatrix<Type> const& transitionMatrix, std::vector<uint_fast64_t> nondeterministicChoiceIndices, storm::storage::SparseMatrix<Type> const& backwardTransitions, storm::storage::BitVector const& initialStates, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, std::shared_ptr<storm::solver::AbstractNondeterministicLinearEquationSolver<Type>> nondeterministicLinearEquationSolver, bool qualitative) {
+                    size_t numberOfStates = phiStates.size();
+                    
                     // We need to identify the states which have to be taken out of the matrix, i.e.
                     // all states that have probability 0 and 1 of satisfying the until-formula.
                     std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01;
                     if (minimize) {
-                        statesWithProbability01 = storm::utility::graph::performProb01Min(this->getModel(), phiStates, psiStates);
+                        statesWithProbability01 = storm::utility::graph::performProb01Min(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, phiStates, psiStates);
                     } else {
-                        statesWithProbability01 = storm::utility::graph::performProb01Max(this->getModel(), phiStates, psiStates);
+                        statesWithProbability01 = storm::utility::graph::performProb01Max(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, phiStates, psiStates);
                     }
                     storm::storage::BitVector statesWithProbability0 = std::move(statesWithProbability01.first);
                     storm::storage::BitVector statesWithProbability1 = std::move(statesWithProbability01.second);
                     
-                    
                     storm::storage::BitVector maybeStates = ~(statesWithProbability0 | statesWithProbability1);
                     LOG4CPLUS_INFO(logger, "Found " << statesWithProbability0.getNumberOfSetBits() << " 'no' states.");
                     LOG4CPLUS_INFO(logger, "Found " << statesWithProbability1.getNumberOfSetBits() << " 'yes' states.");
                     LOG4CPLUS_INFO(logger, "Found " << maybeStates.getNumberOfSetBits() << " 'maybe' states.");
                     
                     // Create resulting vector.
-                    std::vector<Type> result(this->getModel().getNumberOfStates());
+                    std::vector<Type> result(numberOfStates);
                     
                     // Check whether we need to compute exact probabilities for some states.
-                    if (this->getModel().getInitialStates().isDisjointFrom(maybeStates) || qualitative) {
+                    if (initialStates.isDisjointFrom(maybeStates) || qualitative) {
                         if (qualitative) {
                             LOG4CPLUS_INFO(logger, "The formula was checked qualitatively. No exact probabilities were computed.");
                         } else {
@@ -315,20 +316,20 @@ namespace storm {
                         
                         // 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 = transitionMatrix.getSubmatrix(maybeStates, nondeterministicChoiceIndices);
                         
                         // Get the "new" nondeterministic choice indices for the submatrix.
-                        std::vector<uint_fast64_t> subNondeterministicChoiceIndices = storm::utility::vector::getConstrainedOffsetVector(this->getModel().getNondeterministicChoiceIndices(), maybeStates);
+                        std::vector<uint_fast64_t> subNondeterministicChoiceIndices = storm::utility::vector::getConstrainedOffsetVector(nondeterministicChoiceIndices, maybeStates);
                         
                         // 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 = transitionMatrix.getConstrainedRowSumVector(maybeStates, nondeterministicChoiceIndices, statesWithProbability1, submatrix.getRowCount());
                         
                         // Create vector for results for maybe states.
                         std::vector<Type> x(maybeStates.getNumberOfSetBits());
-                                                
+                        
                         // Solve the corresponding system of equations.
-                        this->nondeterministicLinearEquationSolver->solveEquationSystem(minimize, submatrix, x, b, subNondeterministicChoiceIndices);
+                        nondeterministicLinearEquationSolver->solveEquationSystem(minimize, submatrix, x, b, subNondeterministicChoiceIndices);
                         
                         // Set values of resulting vector according to result.
                         storm::utility::vector::setVectorValues<Type>(result, maybeStates, x);
@@ -339,11 +340,15 @@ namespace storm {
                     storm::utility::vector::setVectorValues<Type>(result, statesWithProbability1, storm::utility::constGetOne<Type>());
                     
                     // Finally, compute a scheduler that achieves the extramal value.
-                    storm::storage::TotalScheduler scheduler = this->computeExtremalScheduler(minimize, false, result);
-                    
+                    storm::storage::TotalScheduler scheduler = computeExtremalScheduler(minimize, transitionMatrix, nondeterministicChoiceIndices, result);
+
                     return std::make_pair(result, scheduler);
                 }
                 
+                std::pair<std::vector<Type>, storm::storage::TotalScheduler> checkUntil(bool minimize, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative) const {
+                    return computeUnboundedUntilProbabilities(minimize, this->getModel().getTransitionMatrix(), this->getModel().getNondeterministicChoiceIndices(), this->getModel().getBackwardTransitions(), this->getModel().getInitialStates(), phiStates, psiStates, this->nondeterministicLinearEquationSolver, qualitative);
+                }
+                
                 /*!
                  * Checks the given formula that is an instantaneous reward formula.
                  *
@@ -453,9 +458,9 @@ namespace storm {
                     storm::storage::BitVector infinityStates;
                     storm::storage::BitVector trueStates(this->getModel().getNumberOfStates(), true);
                     if (minimize) {
-                        infinityStates = std::move(storm::utility::graph::performProb1A(this->getModel(), this->getModel().getBackwardTransitions(), trueStates, targetStates));
+                        infinityStates = std::move(storm::utility::graph::performProb1A(this->getModel().getTransitionMatrix(), this->getModel().getNondeterministicChoiceIndices(), this->getModel().getBackwardTransitions(), trueStates, targetStates));
                     } else {
-                        infinityStates = std::move(storm::utility::graph::performProb1E(this->getModel(), this->getModel().getBackwardTransitions(), trueStates, targetStates));
+                        infinityStates = std::move(storm::utility::graph::performProb1E(this->getModel().getTransitionMatrix(), this->getModel().getNondeterministicChoiceIndices(), this->getModel().getBackwardTransitions(), trueStates, targetStates));
                     }
                     infinityStates.complement();
                     storm::storage::BitVector maybeStates = ~targetStates & ~infinityStates;
@@ -526,8 +531,8 @@ namespace storm {
                     storm::utility::vector::setVectorValues(result, infinityStates, storm::utility::constGetInfinity<Type>());
                     
                     // Finally, compute a scheduler that achieves the extramal value.
-                    storm::storage::TotalScheduler scheduler = this->computeExtremalScheduler(this->minimumOperatorStack.top(), false, result);
-                    
+                    storm::storage::TotalScheduler scheduler = computeExtremalScheduler(this->minimumOperatorStack.top(), this->getModel().getTransitionMatrix(), this->getModel().getNondeterministicChoiceIndices(), result, this->getModel().hasStateRewards() ? &this->getModel().getStateRewardVector() : nullptr, this->getModel().hasTransitionRewards() ? &this->getModel().getTransitionRewardMatrix() : nullptr);
+
                     return std::make_pair(result, scheduler);
                 }
                 
@@ -540,33 +545,33 @@ namespace storm {
                  * @param takenChoices The output vector that is to store the taken choices.
                  * @param nondeterministicChoiceIndices The assignment of states to their nondeterministic choices in the matrix.
                  */
-                storm::storage::TotalScheduler computeExtremalScheduler(bool minimize, bool addRewards, std::vector<Type> const& result) const {
-                    std::vector<Type> temporaryResult(this->getModel().getNondeterministicChoiceIndices().size() - 1);
+                static storm::storage::TotalScheduler computeExtremalScheduler(bool minimize, storm::storage::SparseMatrix<Type> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, std::vector<Type> const& result, std::vector<Type> const* stateRewardVector = nullptr, storm::storage::SparseMatrix<Type> const* transitionRewardMatrix = nullptr) {
+                    std::vector<Type> temporaryResult(nondeterministicChoiceIndices.size() - 1);
                     std::vector<Type> nondeterministicResult(result);
                     storm::solver::GmmxxLinearEquationSolver<Type> solver;
-                    solver.performMatrixVectorMultiplication(this->getModel().getTransitionMatrix(), nondeterministicResult, nullptr, 1);
-                    if (addRewards) {
+                    solver.performMatrixVectorMultiplication(transitionMatrix, nondeterministicResult, nullptr, 1);
+                    if (stateRewardVector != nullptr || transitionRewardMatrix != nullptr) {
                         std::vector<Type> totalRewardVector;
-                        if (this->getModel().hasTransitionRewards()) {
-                            std::vector<Type> totalRewardVector = this->getModel().getTransitionMatrix().getPointwiseProductRowSumVector(this->getModel().getTransitionRewardMatrix());
-                            if (this->getModel().hasStateRewards()) {
+                        if (transitionRewardMatrix != nullptr) {
+                            totalRewardVector = transitionMatrix.getPointwiseProductRowSumVector(*transitionRewardMatrix);
+                            if (stateRewardVector != nullptr) {
                                 std::vector<Type> stateRewards(totalRewardVector.size());
-                                storm::utility::vector::selectVectorValuesRepeatedly(stateRewards, storm::storage::BitVector(this->getModel().getStateRewardVector().size(), true), this->getModel().getNondeterministicChoiceIndices(), this->getModel().getStateRewardVector());
+                                storm::utility::vector::selectVectorValuesRepeatedly(stateRewards, storm::storage::BitVector(stateRewardVector->size(), true), nondeterministicChoiceIndices, *stateRewardVector);
                                 storm::utility::vector::addVectorsInPlace(totalRewardVector, stateRewards);
                             }
                         } else {
                             totalRewardVector.resize(nondeterministicResult.size());
-                            storm::utility::vector::selectVectorValuesRepeatedly(totalRewardVector, storm::storage::BitVector(this->getModel().getStateRewardVector().size(), true), this->getModel().getNondeterministicChoiceIndices(), this->getModel().getStateRewardVector());
+                            storm::utility::vector::selectVectorValuesRepeatedly(totalRewardVector, storm::storage::BitVector(stateRewardVector->size(), true), nondeterministicChoiceIndices, *stateRewardVector);
                         }
                         storm::utility::vector::addVectorsInPlace(nondeterministicResult, totalRewardVector);
                     }
                     
-                    std::vector<uint_fast64_t> choices(this->getModel().getNumberOfStates());
+                    std::vector<uint_fast64_t> choices(result.size());
                     
                     if (minimize) {
-                        storm::utility::vector::reduceVectorMin(nondeterministicResult, temporaryResult, this->getModel().getNondeterministicChoiceIndices(), &choices);
+                        storm::utility::vector::reduceVectorMin(nondeterministicResult, temporaryResult, nondeterministicChoiceIndices, &choices);
                     } else {
-                        storm::utility::vector::reduceVectorMax(nondeterministicResult, temporaryResult, this->getModel().getNondeterministicChoiceIndices(), &choices);
+                        storm::utility::vector::reduceVectorMax(nondeterministicResult, temporaryResult, nondeterministicChoiceIndices, &choices);
                     }
 
                     return storm::storage::TotalScheduler(choices);
diff --git a/src/models/AtomicPropositionsLabeling.h b/src/models/AtomicPropositionsLabeling.h
index 486997829..2de73bd1d 100644
--- a/src/models/AtomicPropositionsLabeling.h
+++ b/src/models/AtomicPropositionsLabeling.h
@@ -289,7 +289,7 @@ public:
 	 */
 	void addState() {
 		for(uint_fast64_t i = 0; i < apsCurrent; i++) {
-			singleLabelings[i].resize(singleLabelings[i].getSize() + 1);
+			singleLabelings[i].resize(singleLabelings[i].size() + 1);
 		}
 		stateCount++;
 	}
diff --git a/src/models/Dtmc.h b/src/models/Dtmc.h
index b573f848d..563eee78e 100644
--- a/src/models/Dtmc.h
+++ b/src/models/Dtmc.h
@@ -147,13 +147,13 @@ public:
 		}
 
 		// Does the vector have the right size?
-		if(subSysStates.getSize() != this->getNumberOfStates()) {
+		if(subSysStates.size() != this->getNumberOfStates()) {
 			LOG4CPLUS_INFO(logger, "BitVector has wrong size. Resizing it...");
 			subSysStates.resize(this->getNumberOfStates());
 		}
 
 		// Test if it is a proper subsystem of this Dtmc, i.e. if there is at least one state to be left out.
-		if(subSysStates.getNumberOfSetBits() == subSysStates.getSize()) {
+		if(subSysStates.getNumberOfSetBits() == subSysStates.size()) {
 			LOG4CPLUS_INFO(logger, "All states are kept. This is no proper subsystem.");
 			return storm::models::Dtmc<T>(*this);
 		}
diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h
index ea6c0a3b4..4b3cc1efb 100644
--- a/src/storage/BitVector.h
+++ b/src/storage/BitVector.h
@@ -675,8 +675,8 @@ public:
 	/*!
 	 * Retrieves the number of bits this bit vector can store.
 	 */
-	uint_fast64_t getSize() const {
-		return bitCount;
+	size_t size() const {
+		return static_cast<size_t>(bitCount);
 	}
 
 	/*!
diff --git a/src/utility/graph.h b/src/utility/graph.h
index c9c58a86a..4a7a6f067 100644
--- a/src/utility/graph.h
+++ b/src/utility/graph.h
@@ -177,9 +177,11 @@ namespace storm {
              * @return A bit vector that represents all states with probability 0.
              */
             template <typename T>
-            storm::storage::BitVector performProbGreater0E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0) {
+            storm::storage::BitVector performProbGreater0E(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0) {
+                size_t numberOfStates = phiStates.size();
+                
                 // Prepare resulting bit vector.
-                storm::storage::BitVector statesWithProbabilityGreater0(model.getNumberOfStates());
+                storm::storage::BitVector statesWithProbabilityGreater0(numberOfStates);
                 
                 // Add all psi states as the already satisfy the condition.
                 statesWithProbabilityGreater0 |= psiStates;
@@ -191,9 +193,9 @@ namespace storm {
                 std::vector<uint_fast64_t> stepStack;
                 std::vector<uint_fast64_t> remainingSteps;
                 if (useStepBound) {
-                    stepStack.reserve(model.getNumberOfStates());
+                    stepStack.reserve(numberOfStates);
                     stepStack.insert(stepStack.begin(), psiStates.getNumberOfSetBits(), maximalSteps);
-                    remainingSteps.resize(model.getNumberOfStates());
+                    remainingSteps.resize(numberOfStates);
                     for (auto state : psiStates) {
                         remainingSteps[state] = maximalSteps;
                     }
@@ -230,6 +232,13 @@ namespace storm {
                 return statesWithProbabilityGreater0;
             }
             
+            template <typename T>
+            storm::storage::BitVector performProb0A(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                storm::storage::BitVector statesWithProbability0 = performProbGreater0E(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, phiStates, psiStates);
+                statesWithProbability0.complement();
+                return statesWithProbability0;
+            }
+            
             /*!
              * Computes the sets of states that have probability 0 of satisfying phi until psi under all
              * possible resolutions of non-determinism in a non-deterministic model. Stated differently,
@@ -246,11 +255,9 @@ namespace storm {
              */
             template <typename T>
             storm::storage::BitVector performProb0A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
-                storm::storage::BitVector statesWithProbability0 = performProbGreater0E(model, backwardTransitions, phiStates, psiStates);
-                statesWithProbability0.complement();
-                return statesWithProbability0;
+                return performProb0A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), backwardTransitions, phiStates, psiStates);
             }
-            
+
             /*!
              * Computes the sets of states that have probability 1 of satisfying phi until psi under at least
              * one possible resolution of non-determinism in a non-deterministic model. Stated differently,
@@ -264,15 +271,13 @@ namespace storm {
              * @return A bit vector that represents all states with probability 1.
              */
             template <typename T>
-            storm::storage::BitVector performProb1E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
-                // Get some temporaries for convenience.
-                storm::storage::SparseMatrix<T> const& transitionMatrix = model.getTransitionMatrix();
-                std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
+            storm::storage::BitVector performProb1E(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                size_t numberOfStates = phiStates.size();
                 
                 // Initialize the environment for the iterative algorithm.
-                storm::storage::BitVector currentStates(model.getNumberOfStates(), true);
+                storm::storage::BitVector currentStates(numberOfStates, true);
                 std::vector<uint_fast64_t> stack;
-                stack.reserve(model.getNumberOfStates());
+                stack.reserve(numberOfStates);
                 
                 // Perform the loop as long as the set of states gets larger.
                 bool done = false;
@@ -323,6 +328,15 @@ namespace storm {
                 return currentStates;
             }
             
+            template <typename T>
+            std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01Max(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                std::pair<storm::storage::BitVector, storm::storage::BitVector> result;
+                
+                result.first = performProb0A(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, phiStates, psiStates);
+                result.second = performProb1E(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, phiStates, psiStates);
+                return result;
+            }
+
             /*!
              * Computes the sets of states that have probability 0 or 1, respectively, of satisfying phi
              * until psi in a non-deterministic model in which all non-deterministic choices are resolved
@@ -335,14 +349,7 @@ namespace storm {
              */
             template <typename T>
             std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01Max(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
-                std::pair<storm::storage::BitVector, storm::storage::BitVector> result;
-                
-                // Get the backwards transition relation from the model to ease the search.
-                storm::storage::SparseMatrix<T> backwardTransitions = model.getBackwardTransitions();
-                
-                result.first = performProb0A(model, backwardTransitions, phiStates, psiStates);
-                result.second = performProb1E(model, backwardTransitions, phiStates, psiStates);
-                return result;
+                return performProb01Max(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), model.getBackwardTransitions(), phiStates, psiStates);
             }
             
             /*!
@@ -360,13 +367,11 @@ namespace storm {
              * @return A bit vector that represents all states with probability 0.
              */
             template <typename T>
-            storm::storage::BitVector performProbGreater0A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0) {
-                // Prepare resulting bit vector.
-                storm::storage::BitVector statesWithProbabilityGreater0(model.getNumberOfStates());
+            storm::storage::BitVector performProbGreater0A(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool useStepBound = false, uint_fast64_t maximalSteps = 0) {
+                size_t numberOfStates = phiStates.size();
                 
-                // Get some temporaries for convenience.
-                storm::storage::SparseMatrix<T> const& transitionMatrix = model.getTransitionMatrix();
-                std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
+                // Prepare resulting bit vector.
+                storm::storage::BitVector statesWithProbabilityGreater0(numberOfStates);
                 
                 // Add all psi states as the already satisfy the condition.
                 statesWithProbabilityGreater0 |= psiStates;
@@ -378,9 +383,9 @@ namespace storm {
                 std::vector<uint_fast64_t> stepStack;
                 std::vector<uint_fast64_t> remainingSteps;
                 if (useStepBound) {
-                    stepStack.reserve(model.getNumberOfStates());
+                    stepStack.reserve(numberOfStates);
                     stepStack.insert(stepStack.begin(), psiStates.getNumberOfSetBits(), maximalSteps);
-                    remainingSteps.resize(model.getNumberOfStates());
+                    remainingSteps.resize(numberOfStates);
                     for (auto state : psiStates) {
                         remainingSteps[state] = maximalSteps;
                     }
@@ -452,7 +457,14 @@ namespace storm {
              */
             template <typename T>
             storm::storage::BitVector performProb0E(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
-                storm::storage::BitVector statesWithProbability0 = performProbGreater0A(model, backwardTransitions, phiStates, psiStates);
+                storm::storage::BitVector statesWithProbability0 = performProbGreater0A(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), backwardTransitions, phiStates, psiStates);
+                statesWithProbability0.complement();
+                return statesWithProbability0;
+            }
+
+            template <typename T>
+            storm::storage::BitVector performProb0E(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices,  storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                storm::storage::BitVector statesWithProbability0 = performProbGreater0A(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, phiStates, psiStates);
                 statesWithProbability0.complement();
                 return statesWithProbability0;
             }
@@ -470,15 +482,13 @@ namespace storm {
              * @return A bit vector that represents all states with probability 0.
              */
             template <typename T>
-            storm::storage::BitVector performProb1A(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
-                // Get some temporaries for convenience.
-                storm::storage::SparseMatrix<T> const& transitionMatrix = model.getTransitionMatrix();
-                std::vector<uint_fast64_t> const& nondeterministicChoiceIndices = model.getNondeterministicChoiceIndices();
+            storm::storage::BitVector performProb1A( storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                size_t numberOfStates = phiStates.size();
                 
                 // Initialize the environment for the iterative algorithm.
-                storm::storage::BitVector currentStates(model.getNumberOfStates(), true);
+                storm::storage::BitVector currentStates(numberOfStates, true);
                 std::vector<uint_fast64_t> stack;
-                stack.reserve(model.getNumberOfStates());
+                stack.reserve(numberOfStates);
                 
                 // Perform the loop as long as the set of states gets smaller.
                 bool done = false;
@@ -528,6 +538,15 @@ namespace storm {
                 return currentStates;
             }
             
+            template <typename T>
+            std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01Min(storm::storage::SparseMatrix<T> const& transitionMatrix, std::vector<uint_fast64_t> const& nondeterministicChoiceIndices, storm::storage::SparseMatrix<T> const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
+                std::pair<storm::storage::BitVector, storm::storage::BitVector> result;
+                
+                result.first = performProb0E(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, phiStates, psiStates);
+                result.second = performProb1A(transitionMatrix, nondeterministicChoiceIndices, backwardTransitions, phiStates, psiStates);
+                return result;
+            }
+
             /*!
              * Computes the sets of states that have probability 0 or 1, respectively, of satisfying phi
              * until psi in a non-deterministic model in which all non-deterministic choices are resolved
@@ -540,14 +559,7 @@ namespace storm {
              */
             template <typename T>
             std::pair<storm::storage::BitVector, storm::storage::BitVector> performProb01Min(storm::models::AbstractNondeterministicModel<T> const& model, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates) {
-                std::pair<storm::storage::BitVector, storm::storage::BitVector> result;
-                
-                // Get the backwards transition relation from the model to ease the search.
-                storm::storage::SparseMatrix<T> backwardTransitions = model.getBackwardTransitions();
-                
-                result.first = performProb0E(model, backwardTransitions, phiStates, psiStates);
-                result.second = performProb1A(model, backwardTransitions, phiStates, psiStates);
-                return result;
+                return performProb01Min(model.getTransitionMatrix(), model.getNondeterministicChoiceIndices(), model.getBackwardTransitions(), phiStates, psiStates);
             }
             
             /*!